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ABSTRACT 


Four  topics  are  discussed  in  this  report.  The  first  topic 
is  the  detection  of  underground  nuclear  explosions  in  the 
presence  of  large  natural  events.  Previous  studies  of  this 
problem  have  been  extended  and  partially  confirmed  by  applying 
calculations  to  additional  data  and  by  considering  the  effect  of 
strengthening  certain  assumptions.  The  use  of  a  LASA  for  esti¬ 
mating  epicenters  by  beamforming  and  other  methods  'is  the  second 
topic  considered  in  this  report.  Recent  studies  of  this  topic 
have  employed  computer-aided  simulations  with  controlled  amounts 
of  additive  noise  and  travel-time  anomalies  to  examine  the  rela¬ 
tive  performance  of  different  algorithms  for  epicenter  estimation. 
Simulations  performed  to  date  have  indicated  that  a  beamsplitting 
technique  using  a  DIMUS  array  (i.e.,  only  one  bit  per  sample  on 
each  seismogram)  leads  to  nearly  the  same  level  of  performance 
available  from  the  full  analog  array,  and  both  beamsplitting 
techniques  perform  significantly  better  than  techniques  based 
on  individual  time  picks  on  each  channel,  especially  in  the  case 
of  poor  signal-to-noise  ratios.  The  third  topic  considered  in 
this  report  is  an  analysis  of  detection  and  false  alarm  proba¬ 
bilities  for  a  DIMUS  array.  Taking  advantage  of  a  strong  analogy 
between  this  problem  and  a  standard  problem  in  coding  theory, 
bounds  and  asymptotic  expressions  for  these  probabilities  are 
described,  and  the  signal-to-noise  ratio  loss  of  the  DIMUS  array 
over  the  corresponding  analog  array  is  analyzed  as  a  function  of 
several  parameters.  Finally,  the  previously  developed  automatic 
pP  test  has  been  re-examined  with  a  view  toward  evaluating  its 
performance  and  simplifying  its  implementation. 
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SECTION  I 


INTRODUCTION 

This  is  the  Third  Quarterly  Technical  Report  on  Con¬ 
tract  F19628-67-C-0370 .  The  report  summarizes  progress  on 
four  topics:  the  detection  of  underground  nuclear  explosions 
in  the  presence  of  large  natural  events,  beamforming  and  other 
techniques  for  epicenter  locations,  detection  and  false  alarm 
probabilities  for  DIMUS  arrays,  and  a  re-examination  of  an 
automatic  pP  test. 

Section  II  discusses  the  problem  of  detecting  under¬ 
ground  nuclear  explosions  in  the  presence  of  large  natural 
events.  The  work  described  in  this  section  is  an  extension  of 
previously  reported  work.  In  our  Second  Quarterly  Report,  data 
were  presented  that  indicated  the  sensitivity  of  continental- 
size  arrays  in  detecting  covert  underground  nuclear  tests. 

These  calculations  were  based  on  simulations  using  actual  seismic 
data  from  a  single  large  earthquake.  The  calculations  have  now 
been  repeated  using  data  from  another  earthquake.  These  results 
are  presented  in  Section  II  and  indicate  that  the  previously 
reported  statistics  may  be  regarded  as  representative.  In  our 
First  Quarterly  Report ,  we  suggested  that  a  network  of  single 
seismometers  would  be  useful  in  monitoring  for  a  covert  under¬ 
ground  nuclear  test.  In  that  report  it  was  assumed  that  if  any 
of  the  single  seismometers  received  the  signal  from  an  under¬ 
ground  test  prior  to  the  signal  from  the  earthquake,  the  test 
would  be  detected.  Based  on  this  assumption, minimum  time 
delays  were  calculated  as  a  function  of  the  position  of  the 
test  site  relative  to  the  earthquake  epicenter.  In  this  section, 
we  consider  the  effect  of  strengthening  this  assumption  to  require 
that  at  least  two  or  at  least  three  of  the  single  seismometers 
receive  the  test  signal  before  the  earthquake  signal.  The  con¬ 
sequences  of  strengthening  this  assumption  are  examined  for  two 
collections  of  stations.  Based  on  these  calculations,  it  would 
appear  that  some  upgrading  of  the  quality  of  existing  stations 
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would  be  necessary  to  provide  a  satisfactory  network  of  single 
seismometers  for  this  purpose. 

Section  III  reports  the  results  of  continued  studies 
on  beamforming  and  other  methods  of  epicenter  location  with  a 
LASA.  Our  Second  Quarterly  Report  described  some  theoretical 
analyses  and  simulation  studies  of  this  problem.  The  simulation 
studies  had  been  continued,  with  some  modifications,  during  this 
quarter  and  some  additional  theoretical  analyses  have  been  per¬ 
formed.  It  has  been  found  that  the  automatic  detector  previously 
employed  triggers  earlier  than  desired  on  very  large  events,  and, 
as  a  consequence,  the  location  errors  presented  in  the  last 
report  were  occasionally  larger  than  expected.  The  automatic 
detector  has  been  modified,  and  numerous  simulations  have  been 
performed  with  controlled  amounts  of  additive  noise  and  travel¬ 
time  anomalies.  Pour  different  algorithms  for  estimating  epi¬ 
centers  were  used  in  these  simulations.  Of  these,  the  beam¬ 
splitting  techniques  performed  significantly  better  than  techniques 
based  on  individual  time  picks  on  each  seismogram  whenever  the 
signal-to-noise  ratio  was  small.  Of  particular  interest  in  this 
study  is  the  performance  level  that  has  been  achieved  by  the 
DIMUS  beamsplitting  algorithm,  which  uses  only  one  bit  per  sample 
of  the  individual  seismograms.  In  the  case  of  low  signal-to- 
noise  ratios,  this  algorithm  has  yielded  a  level  of  performance 
almost  equal  to  that  of  the  analog  beamsplitting  techniques, 
and  significantly  better  than  that  available  by  plane-wave 
fitting  after  individual  arrival  time  measurements.  When  the 
signal-to-noise  ratio  is  very  high,  of  course,  there  is  relatively 
little  difference  between  the  different  algorithms.  These  results 
must  be  regarded  as  preliminary,  since  they  are  based  on  only  one 
collection  of  seismic  waveforms  plus  varying  amounts  of  additive 
noise  and  time  anomalies.  However,  it  is  expected  that  similar 
results  will  obtain  when  additional  data  are  considered. 

In  Section  IV,  the  detection  and  false  alarm  probabilities 
for  a  DIMUS  array  are  analyzed  for  a  variety  of  parameter  values. 

A  very  strong  analogy  between  this  problem  and  a  familiar  problem 
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in  coding  theory  is  employed  in  order  to  obtain  fairly  tight 
bounds  on  detection  and  false  alarm  probabilities  as  well  as 
asymptotic  expressions  for  their  behavior.  These  bounds  and 
the  asymptotic  expressions  provide  relatively  easy  means  of 
obtaining  good  estimates  of  these  probabilities  without  having 
to  evaluate  cumbersome  sums  of  binomial  coefficients.  These 
techniques  are  illustrated  in  several  ways  for  various  para¬ 
meter  values, and  the  performance  levels  available  from  a  DIMUS 
array  are  contrasted  with  those  available  from  the  corres¬ 
ponding  analog  array.  It  is  observed  that  in  a  limited  range 
of  parameter  values  the  signal-to-noise  ratio  loss  of  the  DIMUS 
array  compared  to  the  analog  array  is  2  dB ;  in  general,  it 
is  somewhat  larger.  These  results  are  discussed  and  illustrated 
with  specific  examples  chosen  to  be  relevant  to  present  and 
future  seismic  arrays. 

Section  V  presents  a  summary  of  a  re-examination  of 
the  automatic  pP  test,  which  has  been  described  in  earlier 
reports.  This  study  was  undertaken  to  understand  more  clearly 
under  what  circumstances  this  test  yields  results  that  would 
otherwise  be  unavailable,  to  understand  what  irreducible  limita¬ 
tions  on  the  test  performance  do  exist,  and  to  consider  what 
simplifications  in  the  implementation  of  the  test  might  be 
possible  without  a  serious  degradation  in  its  performance.  In 
order  to  consider  these  topics,  several  specific  events  have  been 
reprocessed  and  the  results  are  presented  in  this  section.  In 
addition,  previous  calculations  have  been  reviewed  and  certain 
statistics  abstracted  from  these.  Based  on  these  studies  it 
is  concluded  that  the  use  of  very  large  arrays  for  enhancing 
the  pP  phase  does  yield  results  that  are  sometimes  not  available 
by  individual  examination  of  the  seismograms.  The  adjustment 
for  P-pP  delay  times  appears  to  be  a  useful  procedure  in  com¬ 
bining  seismograms  from  these  very  large  arrays.  In  the  case 
of  moderate  depth  earthquakes  (10-40  km)  it  now  appears  that  the 
fine-grained  calculation  of  the  average  paired  correlation 
coefficient  may  be  unnecessary  and  that  the  same  information 
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could  be  obtained  by  calculating  this  statistic  on  a  10-second 
interval,  thereby  reducing  the  number  of  computations  necessary 
to  carry  out  the  test.  Whether  or  not  this  conclusion  applies 
when  the  test  is  used  for  depths  greater  than  40  km  is  not  yet 
clear.  Finally,  it  is  observed  that  apparent  reflections  resulting 
from  velocity  discontinuities  in  the  earth's  crust  may  produce 
additional  signals  that  cannot  easily  be  distinguished  from  the 
legitimate  pP  and  sP  phases.  It  may  well  be  that  the  only  way 
to  avoid  ambiguities  as  a  result  of  these  reflections  is  to 
acquire  adequate  knowledge  of  the  geology  local  to  the  epicenter 
region . 


SECTION  II 


DETECTION  OF  UNDERGROUND  NUCLEAR  EXPLOSIONS  IN 
THE  PRESENCE  OF  LARGE  NATURAL  EVENTS 

The  problem  of  detecting  an  underground  nuclear  explosion 
that  occurs  shortly  after  an  earthquake  and  in  the  same  geographical 
vicinity  has  been  considered  in  earlier  reports  [1,2].  A  detection 
system  consisting  of  a  continental-size  array  plus  a  set  of  single 
stations  encircling  the  earthquake  epicenter  at  close  ranges  has 
been  considered,  and  the  properties  of  this  hypothetical  system 
have  been  studied.  The  single  stations  impose  a  minimum  delay 
time  between  earthquake  and  test, and  the  continental-size  array 
imposes  a  maximum.  If  either  constraint  is  violated,  the  test 
risks  detection.  If  the  minimum  delay  time  is  violated,  the  signal 
from  the  test  will  reach  at  least  one  of  the  single  seismometers 
before  that  of  the  earthquake.  If  the  maximum  time  delay  is  vio¬ 
lated,  the  earthquake  signals  will  have  decayed  sufficiently  that 
the  continental-size  array  will  be  able  to  detect  the  test.  Results 
of  additional  studies  of  both  components  of  this  detection  system 
are  reported  in  this  section.  Since  this  work  is  a  continuation  of 
earlier  studies  that  are  reported  in  our  First  and  Second  Quarterly 
Technical  Reports,  several  references  to  these  reports  will  appear 
in  the  following  discussion. 

In  the  Second  Quarterly  Technical  Report  [2],  the  depen¬ 
dence  of  the  detection  sensitivity  of  the  array  on  the  time  after 
the  earthquake  and  the  location  of  the  covert  test  relative  to  the 
earthquake  epicenter  was  studied  using  seismic  data  from  a  single 
large  earthquake  recorded  at  fifteen  stations  located  in  the  con¬ 
tinental  United  States.  In  order  to  be  sure  that  these  results 
are  representative,  similar  calculations  have  been  repeated  using 
data  from  another  large  earthquake  and  the  results  of  these  calcu¬ 
lations  are  summarized  below.  Based  on  these  additional  calcu¬ 
lations,  it  is  concluded  that  the  previously  reported  results  may 
be  regarded  as  representative. 
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The  second  component  of  the  proposed  detection  system  is 
a  network  of  single  stations  encircling  the  earthquake  epicenter 
at  close  ranges  (perhaps  40°).  Assuming  that  a  covert  underground 
shot  would  be  detected  if  the  signal  from  the  explosion  reached 
any  of  the  single  seismometers  before  the  earthquake  signal,  a 
minimum  delay  time  between  earthquake  and  shot  may  be  calculated 
for  possible  test  sites  in  the  vicinity  of  the  earthquake.  Results 
of  these  calculations  were  summarized  in  the  First  Quarterly 
Technical  Report  [l]  in  the  form  of  minimum  earthquake-test  delay 
time  contours  on  a  map  of  the  earthquake  epicenter  region. 

Recently,  assumptions  more  realistic  than  the  original  one  that  a 
single  station  receiving  the  explosion  signal  before  the  earth¬ 
quake  signal  would  be  sufficient  to  indicate  the  possibility  of 
a  covert  underground  test  have  been  considered.  In  particular, 
the  requirement  that  two,  or  even  three,  stations  receive  the 
signal  from  the  explosion  before  the  signal  from  the  earthquake 
has  been  considered.  With  the  small  number  of  stations  ori¬ 
ginally  considered,  this  stronger  requirement  significantly 
enlarges  the  geographical  area  covered  by  a  given  minimum  delay 
contour.  If,  however,  a  larger  number  of  stations  is  assumed 
these  contours  can  be  reduced  to  a  more  reasonable  size.  This 
point  will  be  illustrated  in  the  discussion  below  by  presenting 
calculations  based  on  existing  sites  of  present  world-wide  stan¬ 
dard  seimograph  stations. 

2.1  SENSITIVITY  OF  THE  CONTINENTAL-SIZE  ARRAY 

Two  detection  statistics  that  are  based  on  the  records 
available  from  a  continental-size  array  were  proposed  and  their 
properties  were  studied  in  the  Second  Quarterly.  These  studies 
were  based  primarily  on  simulations  employing  seismic  data  ob¬ 
tained  at  fifteen  stations  in  the  continental  United  States  from 
a  magnitude-5-5  earthquake  that  occurred  in  the  Fox  Islands  on 
1  September  1964.  The  same  calculations  have  been  performed  on 
seismic  data  obtained  from  fourteen  stations  in  the  continental 
United  States  from  a  magnitude-6.0  earthquake  that  occurred  in 
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the  Aleutians  on  17  March  1965.  Portions  of  the  results  of 
these  calculations  will  be  summarized  here  with  references  to 
the  relevant  figures  of  the  Second  Quarterly. 

The  two  detection  statistics  under  consideration  may 
be  defined  as  follows.  The  energy  ratio  detection  statistic, 
c ,  is  defined  as  the  average  power  in  the  array  output  in  a  one- 
second  window  divided  by  that  of  the  preceding  five  seconds. 

The  second  detection  statistic,  n ,  is  obtained  by  multiplying 
e  by  the  average-paired-correlation-coefficient  obtained  for  the 
same  one-second  interval,  using  the  individual  seismic  records. 
For  both  statistics  the  array  output  is  obtained  by  aligning  the 
seismograms  to  "steer"  for  a  given  "grid  point,"  and  calcu¬ 
lating  a  weighted  sum  of  the  aligned  seismograms.  It  has  been 
shown  [2]  that  the  maximal-ratio  weighting  of  the  individual 
seismograms  is  superior  to  an  equal  gain  weighting  in  this 
context,  and,  therefore,  only  this  weighting  is  considered  in 
this  report. 

In  the  Second  Quarterly,  the  maximum  excursions  of  the 
n  statistics  were  tabulated  for  each  grid  point  in  the  vicinity 
of  the  earthquake  epicenter  [see  Figure  2.6  of  Reference  2]. 
Based  on  these  results  it  was  inferred  that  a  detection  thres¬ 
hold  of  2.0  would  be  a  conservative  threshold  in  terms  of  false 
alarms.  The  results  of  the  same  calculations  performed  on  the 
magnitude-6.0  earthquake  are  presented  in  Figure  2.1,  and  it  may 
be  observed  that  the  same  conclusion  still  applies.  As  would 
be  expected,  there  are  several  differences  in  the  numbers  re¬ 
sulting  from  the  two  earthquakes;  however,  the  general  range 
of  values  seems  comparable,  and  does  not  appear  to  be  dependent 
upon  the  magnitude  of  the  earthquake. 

In  the  Second  Quarterly  it  was  demonstrated  that  the 
n  statistic  could  be  expected  to  exceed  the  threshold  of  2.0 
whenever  there  was  an  energy  excursion,  as  measured  by  the  e 
statistic,  in  excess  of  10  dB.  As  a  result  of  this  observation, 
by  simply  examining  the  time  dependence  of  the  e  statistic  for 
each  grid  point,  it  is  possible  to  conclude  how  large  an  under¬ 
ground  test  would  be  necessary,  as  a  function  of  time  and  test- 
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site  position,  in  order  to  allow  detection  by  the  continental- 
size  array.  These  calculations  have  been  performed  with  data 
obtained  from  the  magnitude-6.0  earthquake  and  the  results  are 
comparable  to  those  which  would  be  predicted  by  adding  0.5  mag¬ 
nitude  units  to  the  results  of  the  calculations  based  on  the 
magnitude-5 • 5  earthquake. 

One  way  to  summarize  the  detection  sensitivity  of 
the  continental-size  array  is  in  terms  of  the  maximum  delay 
time  between  earthquake  and  underground  explosion  as  a  func¬ 
tion  of  test  site  that  is  imposed  by  the  continental-size  array. 

Using  data  from  the  magnitude-6.0  earthquake , the  time  intervals  for 
concealing  a  magnitude-5 • 5  and  a  magnitude-5.0  shot  are  pre¬ 
sented  in  Figures  2.2  and  2.3-  (These  figures  are  comparable 
to  those  based  on  the  magnitude-5.5  earthquake  reported  in 
Figures  2.12  and  2.13  of  the  Second  Quarterly.)  As  before,  two 
numbers  are  given  at  each  grid  point.  The  first  is  the  minimum 
delay  time  imposed  by  a  hypothetical  network  of  single  seismo¬ 
meters  located  at  a  range  of  about  40°  from  the  earthquake 
epicenter.  The  second  number  is  the  maximum  delay  time  imposed 
by  the  continental-size  array.  In  the  cases  where  the  minimum 
delay  time  exceeds  the  maximum  delay  time,  an  X  is  indicated 
at  the  grid  point  on  the  figure  to  emphasize  that  an  underground 
nuclear  test  of  the  assumed  magnitude  could  not  be  detonated 
at  that  grid  point  without  risking  detection  by  either  the 
single  station  network  or  the  continental-size  array.  Comparing 
the  time  intervals  available  for  concealment  with  those  obtained 
earlier  from  different  data  indicates  values  that  are  roughly 
comparable,  but  by  no  means  identical.  Considering,  for  example, 
a  nuclear  test  that  is  one-half  magnitude  unit  smaller  than 
the  earthquake  leads  to  maximum  delay  times  based  on  the  magni¬ 
tude-6.0  earthquake  data  that  are  in  most  cases  within  a  factor 
of  two  of  those  obtained  using  the  magnitude-5 • 5  earthquake  data.  (Note 
that  in  both  cases  a  shot  that  is  one-half  magnitude  unit  smaller 
than  the  earthquake  is  assumed,  i.e.,  the  shot  magnitude  is  5-5 
in  the  first  case  and  5-0  in  the  second.)  Similar  observations 
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TIME  INTERVALS  AVAILABLE  FOR  CONCEALING 
(  TEST  MAGNITUDE  5.5  ) 
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FIGURE  2.3 

TIME  INTERVALS  AVAILABLE  FOR  CONCEALING 
(  TEST  MAGNITUDE  5.0  ) 
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apply  to  the  case  where  the  shot  is  one  full  magnitude  unit 
smaller  than  the  earthquake.  Based  on  these  calculations ,  it 
is  concluded  that  the  results  presented  earlier  may  be  regarded 
as  representative  and  the  continental-size  array  together  with 
the  network  of  single  seismometers  seriously  limit  the  time 
interval  during  which  a  nuclear  test  could  be  concealed.  As 
before,  it  should  be  noted  that  the  number  of  computations 
implied  by  this  detection  network  is  very  large,  so  that  there 

still  remains  a  serious  practical  question  as  to  whether  such 

a  network  could  reasonably  be  operated. 

2.2  SINGLE  STATION  DETECTION 

In  the  First  Quarterly  it  was  assumed  that  if  any 
member  of  a  network  of  single  stations  received  the  signal 
from  an  underground  test  before  that  of  the  nearby  earthquake, 
the  test  would  not  escape  detection.  Based  on  this  assumption, 
minimum  time  delays  between  earthquake  and  underground  nuclear 
test  were  calculated  as  a  function  of  test  site  position  rela¬ 
tive  to  the  earthquake  epicenter.  The  calculations  originally 
presented  used  a  set  of  eighteen  WWSSS  stations  and  assumed  an 
earthquake  epicenter  in  Mongolia.  The  resulting  minimum  time 
delay  contours  are  very  roughly  circular  with  a  radius  of  2° , 

3. 5°  and  8°  for  minimum  time  delays  of  15,  30  and  60  seconds, 
respectively  [see  Figure  3.2  of  Reference  1].  If  this  original 
assumption  were  modified  to  require  that  at  least  two  of  the 
stations  receive  the  signal  from  the  shot  before  that  of  the 
earthquake,  then  these  contours  would  grow  appreciably.  Assuming 
the  same  set  of  stations  and  the  same  epicenter,  calculations 
show  that  the  15-second  contour  would  stretch  as  far  as  5  degrees 
in  one  direction  and  the  30-second  contour  as  far  as  10  degrees. 
Therefore,  with  this  particular  collection  of  stations,  to 
require  more  than  one  station  significantly  increases  the  size 
of  the  minimum  time  delay  contours. 

If  the  network  of  single  stations  contained  a  great 
many  stations  that  were  well  distributed  in  azimuth,  an  increase 
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in  the  requirement,  from  at  least  one  station  to  at  least  two 
or  even  three  stations  detecting,  would  not  seriously  affect  the 
contours.  This  observation  is  illustrated  by  the  minimum  delay¬ 
time  contours  given  in  Figure  2.4.  These  contours  are  based  on 
a  very  large  number  ( 11 6 )  of  stations  from  the  existing  WWSSS 
network.  The  15,  30  and  60  second  contours  based  on  the  original 
assumption  that  a  single  station  detects  are  given  in  solid 
lines, and  the  contours  implied  by  the  requirement  that  at  least 
three  stations  detect  are  given  by  dashed  lines  in  Figure  2.4. 

As  would  be  expected  in  this  case  of  a  great  many  stations , 
the  region  covered  by  these  contours  is  moderate  in  size  and  not 
strongly  affected  by  the  more  stringent  requirement  of  three 
stations  detecting.  The  stations  used  in  this  calculation  and 
the  epicenter  under  consideration  are  indicated  in  the  map  of 
Figure  2.5. 

The  collection  of  stations  considered  in  the  above 
calculations  are  unrealistic  in  that  many  may  not  be  capable 
of  detecting  moderate  size  shots  even  without  the  earthquake 
interference.  In  a  sense  the  results  of  assuming  this  complete 
set  of  stations  represent  an  upper  bound  on  performance  which 
could  be  achieved  with  existing  station  sites.  In  order  to 
obtain  a  more  realistic  measure  of  what  level  of  performance 
could  be  obtained  from  existing  WWSSS  stations,  the  following 
exercise  was  performed.  From  the  C&GS  Earthquake  Data  Reports 
seven  earthquakes  with  magnitudes  between  4.7  and  5.0  that 
occurred  in  the  Kuriles-Kamchatka  region  during  1966  were 
selected.  Then  23  stations,  each  of  which  reported  at  least 
four  of  the  seven  earthquakes  under  consideration,  were  chosen. 
These  23  stations  are  indicated  in  the  map  of  Figure  2.6. 

Using  these  23  stations  plus  the  same  epicenter  as  before  (in 
the  Kamchatka-Kuriles  region) ,  contours  of  minimum  delay  time 
between  earthquake  and  test  were  calculated  for  each  of  three 
assumptions:  at  least  one  station  receives  the  explosion 

signal  before  the  earthquake  signal,  at  least  two,  at  least 
three.  The  results  of  these  calculations  are  presented  in 
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FIGURE  2.4 

MAXIMUM  TIME  DELAY  CONTOURS 
116  STATION  NETWORK 
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FIGURE  2.5  MAP  OF  116  STATION  NETWORK 
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FIGURE  2,6  MAP  OF  23  STATION  NETWORK 


Figure  2.7A-C.  From  this  figure  it  may  be  observed,  first  of 
all,  that  these  stations  are  not  well  distributed  in  azimuth  and 
that  therefore  test  sites  bearing  150  degrees  from  the  earth¬ 
quake  epicenter  are  not  seriously  affected  by  this  single  station 
network.  Secondly,  it  may  be  observed  in  going  from  one  to  two 
stations,  and  especially  in  going  on  to  three  stations,  the  con¬ 
tours  increase  significantly  in  size.  For  both  reasons  it  does 
not  appear  that  this  particular  collection  of  23  stations  would 
be  satisfactory  as  a  single  station  detection  network.  In  a 
sense,  this  collection  of  stations  represents  a  lower  bound  on 
the  performance  that  could  be  expected  immediately  with  existing 
WWSSS  stations.  From  this  and  the  previous  figure  it  would 
appear  that  a  satisfactory  network  of  single  stations  could  be 
obtained  using  locations  near  those  of  the  WWSSS  stations  provided 
adequate  (single-station)  detection  sensitivity  could  be  achieved 
at  these  locations. 

2.3  SUMMARY  AND  CONCJUSIONS 

Additional  calculations  relating  to  the  problem  of 
detecting  covert  underground  nuclear  tests  have  been  presented. 
These  calculations  represent  a  continuation  of  earlier  study 
on  this  problem  in  which  a  detection  network  consisting  of  a 
continental-size  array  plus  a  network  of  single  stations  well 
distributed  in  azimuth  about  the  epicenter  was  proposed.  In 
earlier  reports  the  sensitivity  of  the  continental-size  array 
was  evaluated  by  using  actual  seismic  records  from  a  single 
magnitude- 5 • 5  earthquake.  Additional  calculations  based  on 
data  from  a  magnitude- 6 . 0  earthquake  have  been  presented  above 
and  suggest  that  the  original  results  may  be  regarded  as 
representative . 

Previously  reported  properties  of  the  single  station 
network  have  been  reconsidered  in  the  context  of  somewhat  more 
realistic  assumptions.  The  original  assumption  was  that  if 
the  signal  from  the  underground  test  arrived  at  any  of  the 
single  seismometers  before  the  earthquake  signal,  the  underground 
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FIGURE  2.7  A 

23  STATION  TIME  DELAY  CONTOURS 
SINGLE  STATION  CRITERION 
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FIGURE  2.7  B 

23  STATION  TIME  DELAY  CONTOURS 
TWO  STATION  CRITERION 
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FIGURE  2.7C 

23  STATION  TIME  DELAY  CONTOURS 
THREE  STATION  CRITERION 
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test  would  be  detected.  The  effect  of  strengthening  this  assump¬ 
tion  to  require  that  at  least  two  or  at  least  three  seismometers 
receive  the  signal  from  the  explosion  before  that  of  the  earth¬ 
quake  has  been  considered  in  this  report .  It  has  been  demonstrated 
that  with  the  set  of  stations  originally  considered,  or  with  another 
set  of  already  existing  stations  totaling  about  twenty,  strengthening 
this  assumption  significantly  enlarges  the  size  of  the  geographic 
area  contained  within  the  minimum  delay  time  contours.  It  would 
seem  from  a  station  location  point  of  view,  however,  that  it 
would  be  possible  to  operate  a  single  station  network  that  per¬ 
formed  adequately.  This  point  is  illustrated  by  calculations 
based  on  a  rather  large  number  (116)  of  WWSSS  stations  for  which 
the  minimum  delay  time  contours  behave  quite  well.  The  conclusion 
suggested  is  that  at  least  the  locations  of  the  seismometers  are 
adequate  for  the  role  of  the  single  station  network,  but  the 
quality  of  many  of  the  individual  stations  would  need  to  be  sub¬ 
stantially  improved. 

Finally,  it  may  be  concluded  that  the  proposed  detection 
network,  consisting  of  a  continental-size  array  plus  a  network  of 
single  stations  at  relatively  close  ranges  and  well  distributed 
in  azimuth,  would  severely  limit  the  possibilities  for  setting  off 
an  undetected  nuclear  explosion  in  the  vicinity  of  an  earthquake. 

The  time  constraints  implied  by  such  a  network  are  on  the  order 
of  two  or  three  minutes  for  an  explosion  of  one-half  to  one 
magnitude  unit  smaller  size  than  the  earthquake.  This  would 
seem  to  be  a  very  severe  restriction  on  the  potential  test  ban 
violator.  One  problem  that  has  not  been  considered  in  this  dis¬ 
cussion  is  the  practical  problem  of  carrying  out  the  large  number 
of  calculations  necessary  to  implement  this  detection  network. 

It  is  possible  that  the  number  of  calculations  might  be  reduced 

considerably  without  a  serious  degradation  in  sensitivity  of 

the  detection  network,  but  this  question  has  not  been  investigated. 
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SECTION  III 


METHODS  FOR  EPICENTER  LOCATION  WITH  LASA 

This  section  is  devoted  to  the  analytical  and  experimental 
study  of  the  errors  in  epicenter  location  that  result  from  travel¬ 
time  anomalies  and  additive  noise.  The  problem  with  which  we  are 
concerned  is  that  of  estimation.  The  detection  of  the  presence 
of  an  event  plays  a  part  here  only  in  its  role  of  defining  which 
portion  of  the  seismogram  is  to  be  considered  as  containing 
signal  plus  noise.  The  detection  is  performed  by  comparing  a 
DIMUS  sum  to  a  threshold.  The  probabilities  of  detection  and 
false  alarm  for  this  detector  are  discussed  in  Section  IV  of  this 
report . 

The  discussion  begins  with  comments  concerning  the  pro¬ 
perties  of  a  DIMUS  detector  when  the  signal-to-noise  ratio  (SNR) 
is  very  large.  It  continues  with  a  consideration  of  the  problems 
caused  by  differences  in  the  waveshapes  observed  at  the  individual 
seismometers.  Then,  results  of  simulations  are  presented, 
together  with  appropriate  analyses,  which  suggest  that  the  method 
described  as  DIMUS  Energy  Beam  Forming  should  be  further  analyzed. 
The  encouraging  results  obtained  and  presented  herein  are  based  on 
too  small  a  data  sample  to  allow  us  to  draw  definite  conclusions. 

In  weak  noise,  it  will  be  shown  that  the  plane  wave  fit  method 
is  adequate  and  reduces  the  complexity  of  estimation.  On  the 
other  hand,  it  falls  short  of  the  expectations  provided  by  the 
beamforming  algorithms  whenever  the  noise  is  not  weak. 

3.1  BACKGROUND  AND  DEFINITIONS 

The  major  portion  of  this  section  is  based  on  Sections 
III  and  IV  of  our  Second  Quarterly  Technical  Report  [2]. 

The  simulation  of  the  effect  of  anomalies  and  noise  is 
based  on  twenty-one  seismograms  recorded  at  the  center  elements 
of  the  subarrays  at  LASA-Montana  for  the  magnitude-5.8  Kazakh 
event  of  21  November  1965.  The  seismograms  were  aligned  so  as 
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to  simulate  a  true  azimuth  of  due  North.  Then  filtered  noise  and 
anomalies  were  added  prior  to  the  application  of  various  algorithms 
used  for  estimation. 

An  automatic  detector  was  incorporated  in  order  to  define 
the  time  intervals  of  interest  in  a  realistic  manner.  The  inte¬ 
gration  windows  were  two  seconds  in  duration  and  centered  about 
the  detection  time  indicated  by  the  DIMUS  detector.  This  detector 
and  various  location  methods  are  discussed  in  Section  IV  of  the 
Second  Quarterly.  We  repeat  a  perfunctory  description  for  conve¬ 
nience,  with  the  addition  of  DIMUS  Energy  Beamforming. 

3.1.1  Automatic  Detection  of  Arrival  Time 

The  implementation  is  carried  out  with  the  output  of 
the  DIMUS  array.  A  threshold  of  19  is  first  assumed  and  location 
in  time  is  given  by  the  violation  of  the  threshold.  If  no  such 
crossing  occurs,  the  threshold  is  successively  decreased  by  one 
unit  and  the  process  repeated  until  a  successful  detection  is 
reached . 

3.1.2  Plane  Wave  Fit  (PWF) 

First,  we  assume  a  rough  knowledge  of  the  azimuth  (here 
due  North).  Then,  arrival  times  determined  by  the  nearest  zero 
crossings  about  the  time  of  detection  are  taken  and  used  to  fit 
a  least-squares  line  with  the  abscissa  determined  by  the  geo¬ 
graphic  projection  of  the  seismometers  on  a  line  perpendicular 
to  the  assumed  direction  of  travel  of  the  plane  wave. 

3.1.3  Analog  Beamforming 

The  array  is  steered  through  various  angles  and  the 
phased  sums  are  computed.  The  two  methods  of  beamforming  estimate 
the  azimuth  by  considering  the  maximum  value  of  the  peak  of 
the  beam  pattern  in  the  method  denoted  as  Analog  Pattern-Peak 
Signal  (AP/P)  or  the  maximum  value  of  the  energy  in  a  two- 
second  interval  of  the  same  beam  in  the  method  denoted  as 
Analog  Pattern-Energy  (AP/E) . 
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3.1.4  DIMUS  Beamforming 

This  method  differs  from  those  previously  proposed.  The 
DIMUS  beamforming  algorithm  performs  the  same  computations  as  the 
analog  beamforming  method  with  the  replacement  of  the  analog  wave¬ 
forms  by  their  DIMUS  equivalents.  It  is  denoted  by  DP/E. 

3.1.5  Synthetic  and  Actual  Seismograms 

We  intend  to  differentiate  between  the  performance  of 
the  algorithms , under  the  assumptions  of  identical  signals  plus 
independent  noises,  and  their  performance  with  the  actual  waveforms 
recorded..  Synthetic  waveforms  are  obtained  by  constructing 
twenty-one  identical  seismograms,  each  being  the  phased  sum 
of  the  original  ones. 

3.1.6  Assumptions  and  Definitions 

Throughout  the  simulations  reported,  a  constant  phase 
velocity  of  20  km/sec  is  assumed  and  the  signal-to-noise  ratio 
(SNR)  is  defined  as 

■y[Wave  Peak]^ 

Ave .  Noise  Power 

This  definition  is  consistent  with  the  conventional  definition 
of  SNR  as  the  ratio  of  signal  energy  and  average  noise  power 
when  sinusoidal  signals  are  considered.  The  justification  for 
its  choice  is  merely  one  of  convenience. 

The  mean,  p ,  reported  in  the  Tables  represents  the 
average  of  the  errors  made  by  the  various  algorithms.  The 
standard  deviation,  0,  is  that  under  the  a  priori  hypothesis 
of  a  zero  mean,  and  thus  is  computed  as  the  rms  value  of  the 
errors  obtained  in  simulation. 
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3.2 


DIMUS  SUM  DETECTION  WITH  LARGE  SNR 


When  the  DIMUS  sum  is  used  for  detection  of  the  presence 
of  the  signal,  it  is  assumed  that  over  intervals  during  which 
there  is  no  signal  present,  the  number  of  occurrences  of  a  "+1" 
and  a  "-1"  on  each  seismogram  will  average  out  and  not  contribute 
to  any  detectable  congregate.  It  becomes  apparent,  however,  that 
in  the  twenty-one  seismograms  for  the  Kazakh  event  used  in  pre¬ 
viously  reported  simulations  [2],  a  negative  precursor  as  illus¬ 
trated  in  Figure  3 • 1A  is  present  in  all  seismograms.  The  fact 
that  this  precursor  seems  to  be  of  negligible  amplitude  in  the 
analog  waveform  is  obviated  when  the  DIMUS  beam  is  formed.  As 
a  result,  the  main  signal  is  presumed  present  prior  to  its 
principal  motion  and  all  time  windows  or  vernier  measurements 
used  in  the  epicenter  location  algorithms  are  subjected  to  an 
error,  as  indicated  in  Figure  3 • IB . 

This  artifact  is  responsible  for  the  deviations  from 
theoretical  predictions  and  the  nonzero  mean  of  the  estimate 
reported  earlier  for  estimated  azimuth  in  the  presence  of 
anomalies . 

Three  methods  to  relieve  this  difficulty  are  presented 
below  with  their  merits  or  shortcomings.  The  first  two  are 
expedient  measures  and  do  not  purport  to  be  suggestions  for  an 
operational  system. 

3.2.1  Null-Zone  Randomization 

The  analog  waveform  is  first  operated  on  by  a  hard- 
limiter  with  a  null  zone.  The  output  of  this  operation  is  a 
three-level  quantization  representation  of  the  analog  waveform. 

Such  precursors  as  mentioned  above,  or  noise  of  small  amplitude, 
would  therefore  be  represented  as  zeroes  rather  than  by  their 
polarity.  On  the  other  hand,  in  order  to  preserve  the  economy 
of  one  bit  of  information  per  sample  provided  by  DIMUS,  the 
null  zone  values  are  randomized  to  plus  and  minus  one  equiprobably . 
The  degradation  of  this  latter  step  is  minimized  by  the  fact 
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FIGURE  3.1 

TYPICAL  SEISMOGRAM  AND  DIMUS  SUM 
FOR  STRONG  EVENT 
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that  the  random  noise  introduced  is  out-of-band  with  respect 
to  the  signal. 

Simulation  results  in  the  presence  of  anomalies  for  the 
same  event  (Kazakh,  21  November  1965)  are  presented  in  Table  I. 
They  are  seen  to  agree  with  theoretical  predictions  derived  in 
the  Second  Quarterly  and  suggest  that  the  earlier  difficulties 
were  due  to  the  effects  of  "false  triggering"  by  the  automatic 
detector . 

The  assignment  of  the  null  zone  threshold  level  is  a 
potential  problem  in  the  use  of  this  method.  It  may  not  be  set 
too  low  with  respect  to  the  ambient  noise,  for  this  would  offer 
no  help  to  the  problem  in  the  case  of  large  events.  On  the 
other  hand,  if  it  is  set  higher  than  the  ambient  noise,  small 
events  remain  in  the  null  zone  for  their  duration  and  become 
undetectable.  The  answer  may  be  a  posteriori  setting,  i.e., 
a  level  determined  by  the  average  energy  received.  If  a  small 
energy  is  detected,  with  or  without  signal  present,  the  null 
zone  remains  negligible.  If  a  large  energy  is  detected  on  a 
seismogram,  the  relatively  higher  threshold  would  eliminate  the 
precursor.  This  method,  however,  introduces  the  complexity 
of  energy  detection  and  memory  since  DIMUS  can  only  be  provided 
after  the  threshold  is  established. 

3.2.2  Synthetic  Noise 

Another  method  of  alleviating  the  problems  caused  by 
the  small  precursor  consists  of  adding  low  level  out-of-band 
noise  to  the  analog  waveform  prior  to  hardlimiting .  This  would 

have  the  same  effect  on  the  precursor  for  DIMUS  as  the  null  zone 
randomization.  For  strong  events,  the  out-of-band  noise  would  be 
annihilated  by  hardlimiting  when  the  signal  is  present. 

Whenever  the  signal  is  small,  the  out-of-band  noise 
would  be  eliminated  by  the  DIMUS  sum  only  if  the  number  of 
seismograms  were  large.  On  the  other  hand,  if  the  scheme  re¬ 
solves  the  difficulty  of  the  precursor,  then  signals  of  lower 
amplitude  than  the  precursor  become  de  facto  undetectable „ 
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TABLE  I 


RUN  1 


RUN  2 


RUN  3 


RUN  4 


RUN  5 


RUN  6 


RUN  7 


SAMPLE  MEANS  AND  STANDARD  DEVIATIONS  VS.  TIME  ANOMALIES 


20  CASES  -  ACTUAL  SEISMOGRAMS,  ±1  SAMPLE  POINT  ANOMALIES 


Vi 

0 

PWF 

0.0  3° 

0.19° 

AP/E 

0.05 

0.30 

AP/P 

0.07 

0.23 

DP/E 

0.09 

0.23 

20  CASES 

-  SYNTHETIC 

SEISMOGRAMS, 

-1  SAMPLE  POINT  ANOMALIES 

u 

0 

PWF 

0.03° 

0.19° 

AP/E 

0.008 

0.17 

AP/P 

0.10 

0.30 

DP/E 

0.008 

0.27 

11  CASES 

-  SYNTHETIC 

SEISMOGRAMS , 

-1,  SAMPLE  POINT 

y 

0 

PWF 

0.04° 

0.17° 

AP/E 

0.06 

0.11 

AP/P 

0.14 

0.28 

24  CASES 

-  ACTUAL  SEISMOGRAMS,  -1 

SAMPLE  POINT 

V 

O 

PWF 

-0.02° 

0.18° 

AP/E 

0.08 

0.34 

AP/P 

0.10 

0.28 

24  CASES 

-  ACTUAL  SEISMOGRAMS,  ±2 

SAMPLE  POINTS 

u 

0 

PWF 

o 

t>- 

i — 1 

o 

1 

0 . 39° 

24  CASES 

-  ACTUAL  SEISMOGRAMS  -  NC 

i  NULL  ZONE  RANDOMIZATION, 

±1  SAMPLE  POINT 

y 

a 

PWF 

o 

ro 

oo 

o 

o 

C\J 

^r 

o 

11  CASES 

-  ACTUAL  SEISMOGRAMS  -  NC 

i  NULL  ZONE  RANDOMIZATION, 

-2  SAMPLE  POINTS 

u 

a 

PWF 

0 . 49° 

o 

o 

o 
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3.2.3 


Analog  Waveform  Post-Decision 


If  we  assume  that  at  least  one  analog  waveform  is 
available  to  the  processor,  in  addition  to  the  DIMUS  waveforms, 
then  a  comparison  with  the  analog  seismograms,  post-detection, 
would  ascertain  whether  or  not  the  response  of  the  detector  is 
due  to  a  precursor.  Since  the  dilemma  is  present  only  at  high 
SNR,  there  is  no  ambiguity  in  this  post-detection  decision. 
Whenever  the  analog  waveform  indicates  that  the  DIMUS  detector 
"triggered"  on  a  small  precursor,  the  indicated  arrival  time 
of  the  event  could  be  modified  appropriately. 

3.3  WAVE  SHAPE  VARIATIONS  ACROSS  THE  ARRAY 

In  any  of  the  methods  used  in  estimating  the  azimuth, 
plane  wave  fit  using  a  zero  crossing  measure  of  arrival  time, 
plane  wave  fit  using  a  correlation  time  pick,  analog  pattern 
with  energy,  analog  pattern  with  phased  sum  peak,  DIMUS  energy 
pattern,  the  performance  is  degraded  by  the  fact  that  the 
signal  portions  of  the  seismograms  are  not  identical.  Indeed, 
it  is  found  that,  for  Kazakh  events,  the  LASA-Montana  center 
elements  B-3  and  D-l ,  for  instance,  yield  a  disparity  of  one- 
third  of  a  second  in  the  time  of  arrival  obtained  by  zero  crossing 
and  correlation  time  pick,  the  delay  being  observed  in  the  latter 
method . 

At  the  price  of  a  sizable  increase  of  complexity,  pre¬ 
processing  of  the  individual  seismometers  may  be  attempted. 

Let  us  assume,  as  a  perfunctory  examination  of  the  strong  event 
studied  seems  to  suggest,  that  the  seismograms  are  made  up  of 
the  superposition  of  several  delayed  arrivals  of  a  basic  waveform. 
If  we  further  assume  that  the  delays,  following  the  first  occur¬ 
rence,  are  not  coherent,  then  the  analog  phased  sum  represents 
a  first-order  approximation  of  the  basic  waveform.  Attempts  to 
use  this  in  order  to  remove  the  interference  in  the  presence  of 
noise  introduces  errors  which  often  fail  to  remove  the  problem 
and  lead  to  a  performance  of  the  epicenter  location  algorithm 
which  is  no  better  than  the  simpler  decision  to  disregard 
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the  delinquent  seismograms  entirely.  In  order  to  recognize  this 
delinquency,  one  may  either  rely  on  the  fact  that  this  behavior 
seems  to  be  consistent  with  the  area  under  observation  or  merely 
use  two  different  methods  of  time-of-arrival  pick  and  set  a 
criterion  on  the  disparity  between  them  for  the  exclusion  of  the 
seismogram . 

The  effect  of  these  differences  is  shown  in  Table  II 
where  it  may  be  observed  in  five  of  six  cases  that  the  rms  errors 
(a  in  the  table)  were  worse  when  the  actual  waveforms  were  used 
than  when  the  same  waveform  was  used  at  each  element. 


3.4  COMPARISON  OF  PLANE  WAVE  FIT  AND  BEAMFORMING 

Let  us  consider  the  idealized  situation  where 
i)  The  signal,  if  it  exists,  is 


<f>  ( t-t '  ) 


/2WE 

O 


sin2irW ( t-t '  ) 

2ttW  ( t-t '  ) 


(3.1) 


where  t'  is  the  time  of  arrival, 
ii)  The  signal  being  bandlimited  to  (-W,W),  we  assume  the 
noise  to  be  bandlimited  to  the  same  interval  with 
double-sided  power  spectrum  Nq/2.  Let  the  statistics 
of  the  noise  be  Gaussian  with  zero  mean. 


3.4.1  Weak  Noise  Errors  for  Time  Picks 


If  a  plane  wave  fit  is  performed,  using  the  time  picks 
for  arrival  of  the  signals  with  twenty-one  seismometers ,  then 
in  view  of  Equation  (A-8)  of  Reference  [2],  where  now  the  zero 
crossing  occurs  at  t=l/2W,  the  mean  square  error  in  the  time 
picks  is 


var[n( t ) ] 

[s'(iw)]2 


(3.2) 


since  1/2W  is  the  location  of  a  zero  crossing. 


“2 

ez 


WN 


0 


N 


0 


SW^E, 


8W2E, 


(3.3) 


S  w"  ~S 

for  a  zero  crossing  time  pick,  and  from  Equation  (A-7) 
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TABLE  II 


RUN  1  -  10  dB, 

20  CASES 

ACTUAL 

SEISMOGRAMS 

SYNTHETIC 

SEISMOGRAMS 

M 

o 

y 

a 

PWF 

-0.29° 

0.50° 

O 

C\J 

o 

o 

0 .  45° 

AP/E 

1 

O 

o 

VO 

0.-41 

-0.07 

0 . 32 

AP/P 

oo 

o 

o 

1 

0.35 

1 

o 

o 

ro 

0.28 

RUN  2  -  3  dB, 

20  CASES 

ACTUAL 

SEISMOGRAMS 

SYNTHETIC 

SEISMOGRAMS 

y 

a 

y 

a 

PWF 

O 

VO 

o 

1 

0.83 

1 

O 

O 

OO 

0.46 

AP/E 

0.02 

0. 3^ 

0.001 

0.36 

AP/P 

CVJ 

o 

o 

1 

0 . 40 

0.002 

0.24 
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of  Appendix  A 


ec 


12 

2 


N, 


 3 


(4W)2  2ES 

n  N  n 


8-rr 2  W2  ES 


(3.4) 


for  a  correlation  time  pick. 

Thus,  if  we  define  G  as  the  power  gain  in  dB  of  the 
correlation  time  pick  over  the  zero  crossing  time  pick,  we 
obtain  2 

G  =  10  log10  j-  *  5  dB  (3.5) 


3.4.2  Strong  Noise  Errors  for  Time  Picks 

Because  of  its  apparent  superiority,  we  examine  the 
large  noise  behavior  of  the  correlation  time  pick  alone.  It 
seems  intuitively  obvious  that  in  strong  additive  noise,  the 
zero  crossing  time  pick  will  degenerate  far  sooner  than  the 
algorithm  performing  correlation  with  a  replica  of  the  signal 
(or  an  estimate  of  this  replica). 

Consider  Figure  3.2  which  depicts  the  output  of  the 
correlator  used  for  estimating  the  time  of  arrival.  Let  tg  be 
the  true  arrival  time  and  t  the  estimate. 

In  order  to  calculate  an  approximation  to  the  probability 
of  such  an  anomalous  error,  let  us  focus  our  attention  on  a 
finite  set  of  instants  {t^}  where  i  is  an  integer,  defined  as 
follows : 


i  =  -K , . . .  ,+K 

The  reason  for  a  finite  limit  on  i  is  that  it  is  unlikely  that 
the  seismograms  will  be  used  for  all  time  in  order  to  obtain 
an  estimate.  In  other  words,  there  is  a  priori  knowledge  that 
if  an  event  has  been  detected,  the  signal  is  assumed  to  be  in 
a  bounded  time  interval.  Thus  we  consider  only  2K  instants  other 
than  tg.  Typically,  K  will  be  a  small  integer. 
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SIGNAL 


FIGURE  3.2 

ANOMALOUS  SIGNAL  DETECTION 
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For  any  pair  of  instants  in  the  interval,  the  waveforms 
have  the  property 


/  ip(t-t1)*(t-tj)dt  =  6ijEs  (3.7) 


where  6^.  is  the  Kroeneker  delta,  the  set  representing  a  set  of 
sampling  functions.  Thus,  the  probability  that  the  noisy  func¬ 
tion  f(t)  will  peak  at  one  of  the  incorrect  time  instants,  or 
in  its  neighborhood,  in  the  set  {t^}  is  the  probability  of  error 
in  the  communication  of  one  of  (2K+1)  equally  likely  messages 
when  using  equal-energy  orthogonal  signals  [3,  Ch.  VIII].  Let 
us  denote  by  r  the  event 


r  =  [f(t^)>f(tg)  for  at  least  one  t^  in  {t^}] 


(3.8) 


Then  we  have  [3] 

oo 

P  { r }  -  1  -  /  :  e 

-oo  /tjN^ 


(y-/Ep 


N 


[1-Q(-^^)]2K  dY 

/iV2 


(3.9) 


where 


Q(o)  =  /  —  e  e  /2  d3 


a  /2  tv 

Approximating  the  Q  function,  we  obtain 

E0 


P(r  }  * 


2K 


/2u(Es/Nq) 


exp{-  gif2} 


0 


(3.10) 


(3.11) 


Since  the  probability  that  the  noise  causes  f(t)  to 
peak  near  t^  is  the  same  as  the  probability  that  f(t)  peaks 
near  t.  for  t^t.  ,  the  error  will  be  approximately  uniformly 
distributed  on  (-K/2W,  K/2W) .  Thus 


2  <  K4 


-Es/2No 


6W‘ 


/2tt(Es/N0) 


(3.12) 


where  1/2W  represents  the  interval  between  two  adjacent  instants 


and  t 


i+1  ’ 
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For  very  low  SNR,  another  approximation  to  the  Q  func¬ 
tion  yields 

— — .  1| 

e2  -  -  /E„/Nn]  (3.13) 

6U3  2  b  U 

Thus,  if  we  let  be  the  mean  square  error  when  the  bound  on 
the  deviation  is  one  pulse  width  (K=2)  on  either  side  and  the 
correlation  time  pick  is  used. 


2  <  3  1  N0  ,  8  1  VN0 

e  -  — 5-  s—  +  o  ■  •  -  e 

8tt  W  S  3W-*  /2tt(Es/Nq) 


(3-l4) 


3.4.3  Energy  Beamforming 

In  the  Second  Quarterly  it  was  shown  (Appendix  B)  that, 
for  weak  noise,  correlation  time  picks  followed  by  a  plane  wave 
fit  and  energy  beamforming  are  identical.  The  behavior  of  the 
two  methods  diverge  in  stronger  noise.  Apart  from  the  requirement 
of  having  a  replica  of  the  signal  in  order  to  perform  a  corre¬ 
lation  time  pick,  the  superiority  of  beamforming  is  borne  out 
by  the  tightness  of  the  linear  approximation  in  Equations  (A-4) 
and  (B-4)  in  Reference  [2]. 

The  linearization  of  Equation  (A-4)  for  the  time  pick 
is  valid  if  the  estimate  -t  of  the  time  of  arrival  is  close  to 
the  true  value  t,  i.e.,  if  the  variance  (t-t)  is  small.  The 
linearization  of  Equation  (B-4)  is  valid  if  the  variance  of 
the  final  estimate  m  of  m  is  small. 

We  must  bear  in  mind,  however,  that  the  relation  between 
the  two  variances  is  that  of  the  noise  reduction  provided  by  plane 
wave  fit ,  i.e,. 


(m-m)' 


12 


,!.;)2  N (N  — 1 ) 


(3.15) 


Thus,  the  signal-to-noise  ratios  where  either  of  the  analyses 
break  down  bear  the  same  relationship. 


-35- 


If  the  case  for  energy  beamforming  based  on  the  analog 
waveforms  is  well  established,  that  of  energy  beamforming  based 
on  the  DIMUS  waveforms  remains  unsettled.  Table  III  offers  a 
comparison  based  on  simulations  for  several  values  of  SNR.  It 
seems  at  first  glance  that  the  analog  and  DIMUS  beamforming 
perform  with  little  statistical  difference,  at  least  in  the 
limited  amount  of  simulation  presented.  We  now  offer  a  rough 
argument  for  this  similarity. 

3.5  COMPARISON  OP  ANALOG  AND  DIMUS  BEAMFORMING 

We  present  the  argument  for  the  similarity  in  performance 
of  these  methods  of  estimating  the  epicenter  location  at  the 
extremes  of  SNR. 

3-5-1  Strong  Noise  Case 

In  strong  noise,  the  property  of  DIMUS  phased  sums  is 
that  they  yield  a  statistically  consistent  representation  (within 
a  constant)  of  the  analog  phased  sums  of  the  waveforms  on  which 
they  are  based.  Thus  the  performance  of  DIMUS  beamforming  is 
expected  to  be  comparable  to  analog  (within  the  loss  of  ir/2(  2  dB) 
in  SNR  associated  with  DIMUS  sums),  for  the  simple  reason  that 
a  similar  answer  is  expected  in  each  case. 

3-5-2  Weak  Noise  Case 

When  the  SNR  is  large,  the  polarity  of  the  DIMUS 

waveform  obtained  from  a  particular  seismogram,  and  that  obtained 

from  the  noise-free  signal  differ  only  in  the  close  neighborhood 
of  the  zero  crossings  of  the  analog  signal.  This  deviation,  further¬ 
more,  is  a  function  of  the  slope  of  the  signal  at  those  zero- 
crossings  . 

DIMUS  beamforming  then  maximizes , all  else  being  assumed 
equal,  the  average  pair-wise  correlation  of  the  DIMUS  waveforms 
near  the  zero  crossings  of  the  analog  waveforms.  The  number  of 
zero  crossings  involved  depends  on  the  width  of  the  integration 
window  used  in  determining  the  energy  in  the  phased  sum. 
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TABLE  III 

COMPARISON  OF  ANALOG  AND  DIMUS  BEAMFORMING 


RUN  1-18  CASES,  NO  ANOMALIES 


10  dB 
3  dB 
-3  dB 


ANALOG 
y  o 

0.016  0.286 

0.081  0.339 

0.065  0.317 


DIMUS 

y  a 

0.035  0.241 

0.068  0.340 

0.012  0.380 


RUN  2-20  CASES.  NO  NOISE,  ±1  SAMPLE  POINT  ANOMALIES 

ANALOG  DIMUS 

y  a  y  a 

Actual  Seismograms  -O.O58  0.302  -0.093  0.239 


Synthetic  Seismograms  -0.008  0.177  -0.008  0.277 
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This  operation,  however,  is  akin  to  minimizing  the 
deviation  from  a  straight  line  of  the  position  of  these  zero- 
crossings  or  equivalently  the  sign  changes  of  the  DIMUS  wave¬ 
forms,  thereby  maximizing  the  number  of  coherent  points.  The 
similarity  with  plane-wave  fit  using  correlation  time  picks  is 
now  placed  in  evidence. 

The  argument,  however,  requires  a  final  comment.  The 
correlation  time  pick  here  would  be  between  a  DIMUS  representa¬ 
tion  of  the  seismogram  and  the  DIMUS  representation  of  the 
noise-free  signal.  If,  as  assumed,  the  noise  is  small,  the 
identical  correlation  operation  with  analog  waveforms  tests 
their  coherence  of  polarity.  Bringing  this  comment  to  its 
conclusion,  we  expect  now  similarity  of  answer  between  analog 
and  DIMUS  beamforming  in  each  particular  test  case  and,  there¬ 
fore,  a  similarity  in  statistical  performance. 

We  now  must  reiterate  our  description  of  the  argument 
as  being  rough.  We  do,  however,  present  in  Table  IV  the 
detailed  results  of  simulation  and,  leaving  the  question  open, 
suggest  that  some  effort  be  spent  to  further  analyze  the 
performance  of  DIMUS  beamforming  in  view  of  its  benefits  in 
reduction  of  computational  complexity. 

3.6  SUMMARY  AND  CONCLUSIONS 

The  behavior  of  the  error  in  the  presence  of  weak  noise 
with  anomalies,  presented  in  Table  I,  in  conjunction  with 
similar  results  in  Section  IV  of  the  Second  Quarterly  [2]  show 
that  the  accuracy  of  the  location  of  the  time  interval  of  interest 
in  the  seismograms  is  of  major  importance.  The  solution  adopted 
in  the  course  of  the  simulations  was  that  of  null  zone  randomi¬ 
zation.  This  proves  to  be  useful  if  there  is  a  priori  knowledge 
that  the  event  is  of  large  magnitude.  On  the  other  hand,  the 
alternative  of  using  the  method  of  analog  waveform  post-decision 
may  be  better  suited  for  a  system  which  must  work  with  signals 
of  a  wide  range  in  magnitude. 
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TABLE  IV 


ERROR  VALUES  FOR  DIMUS  AND  ANALOG  BEAMFORMING 


10 

dB 

3  dB 

-3  dB 

ANALOG 

DIMUS 

ANALOG 

DIMUS 

ANALOG 

DIMUS 

-0. 325° 

-0 . 300° 

0.300 

0.200 

0 . 300 

0.200 

-0.050 

-0.000 

0.250 

0.250 

0.250 

-0.050 

-0.250 

-0.175 

-0.325 

-0.325 

0.550 

0.325 

0 . 300 

0.250 

-0.300 

-0.300 

O 

o 

on 

o 

1 

-0.300 

0.250 

0 . 300 

0.350 

0.350 

-0.325 

-0.350 

0.250 

0.075 

-0.250 

0 . 300 

0.250 

-0.050 

0.275 

0.200 

0.550 

0.550 

0.350 

0.350 

-0.325 

-0.325 

-0.200 

-0.250 

0.249 

-0.300 

0.550 

0.350 

0.750 

0.125 

-0.175 

-0.175 

0.075 

0.000 

0 . 550 

0.550 

O 

O 

on 

o 

0 . 300 

-0.325 

-0.325 

0.075 

0.125 

0.300 

0 . 300 

0 . 550 

0.550 

-0.325 

0.275 

-0.325 

-0.325 

0.300 

0 . 300 

0.275 

0.350 
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Analysis  and  simulation  of  the  several  methods  of 
locating  the  epicenter  allow  us  to  rank  them  according  to  their 
performance.  In  additive  weak  noise,  although  the  performances 
are  uniformly  good,  we  expect,  from  worst  to  best,  plane  wave  fit 
on  zero  crossing  time  picks,  DIMUS  energy  beamforming,  plane 
wave  fit  on  correlation  time  picks  and  analog  beamforming.  In 
strong  noise,  plane  wave  fit  on  zero  crossing,  as  suggested  by 
the  simulation,  is  the  first  one  to  degenerate.  As  the  SNR  is 
further  decreased,  plane  wave  fit  on  correlation  time  picks  is 
expected  to  fall  behind  the  beamforming  methods. 

The  performance  of  DIMUS  energy  beamforming,  in  compari¬ 
son  with  the  analog  beamforming,  suggests  that  a  large  saving 
in  computational  complexity  may  exist  without  any  appreciable 
loss  in  performance.  The  limited  amount  of  simulation,  together 
with  the  simplified  argument  of  Section  3-5,  without  resolving 
the  question,  do  establish  that  further  study  is  indicated. 
However,  it  must  be  pointed  out  that  the  results  of  simulation 
with  the  Kazakh  event  mentioned  earlier  are  in  accordance  with 
this  conjecture. 
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SECTION  IV 


DETECTION  PROBABILITIES  FOR  DIMUS  ARRAYS 

The  detection  capabilities  of  an  array  of  hardlimited 
waveforms  are  considered  in  this  section.  The  objective  of  this 
study  is  a  clear  understanding  of  the  potential  performance  of 
present  and  future  seismic  arrays  when  operated  as  DIMUS  arrays. 
Since  probabilities  of  detection  and  false  alarms  are  often 
difficult  to  evaluate  exactly,  upper  and  lower  bounds,  as  well 
as  the  asymptotic  form  of  these  probabilities,  will  be  presented. 
These  results  should  provide  a  convenient  means  of  evaluating 
the  detection  performance  of  a  DIMUS  array  as  a  function  of  the 
signal-to-noise  ratio  (SNR)  and  the  number  of  elements. 

This  section  begins  with  a  definition  of  the  seismic 
context  in  which  DIMUS  is  of  potential  interest.  The  DIMUS 
detection  problem  is  then  described  mathematically  and  expres¬ 
sions  are  obtained  for  detection  and  false  alarm  probabilities. 
Borrowing  results  from  coding  theory,  these  expressions  are 
bounded  and  their  asymptotic  behavior  is  examined.  The  corres¬ 
ponding  analog  array  is  then  compared  to  the  DIMUS  array  and 
it  is  shown  that  the  SNR  loss  is  at  least  2  dB,  and  may  be 
considerably  larger.  These  results  are  illustrated  in  several 
ways,  including  standard  vs.  SNR  plots  for  various  choices 
of  false  alarm  rate  and  number  of  elements.  Finally,  a  brief 
summary  of  the  relationship  between  the  DIMUS  detection  problem 
and  coding  theory  is  presented. 

M.l  DEFINITIONS  AND  ASSUMPTIONS 

The  idealized  seismic  problem  of  interest  is  confined 

to  a  single  sample  point.  It  is  assumed  that  prior  to  hard- 

limiting  each  of  N  elements  has  a  "signal"  component  of  amplitude 

s  and  a  noise  component,  n^,  that  is  a  zero  mean  Gaussian  random 

variable.  It  is  further  assumed  that  the  n.  are  independent  between 

2  ^ 

elements  and  have  the  same  variance,  o  .  The  DIMUS  array  output  is 
then  given  by 
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3 


N 

l  sgn( s+n. ) 
i-1 

whereas  the  corresponding  "analog"  array  output  would  be  given  by 
N 

l  (s+n. ) . 
i-1 

The  detection  problem  of  interest  is  that  of  deciding  whether 
or  not  the  "signal"  is  present  (i.e.,  s/0).  It  is  easy  to  show 
under  a  variety  of  assumptions  that  the  optimum  decision  rule  is 
simply  to  compare  the  array  output  (DIMUS  or  analog)  to  a  threshold, 
which  is  determined  by  specifying  a  false  alarm  rate,  by  assuming 
a  cost  structure,  or  in  some  other  way. 

We  are  interested  in  evaluating  the  detection  and  false 

alarm  probabilities  as  a  function  of  the  threshold  value,  the  SNR 

2  2 

(e.g.,  s  /a  ) ,  and  the  number  of  elements,  N.  In  addition  we  wish 
to  compare  the  performance  of  the  DIMUS  and  analog  arrays.  It  turns 
out  that  there  is  a  strong  analogy  between  this  problem  and  a 
familiar  one  in  coding  theory.  This  analogy  yields  many  results 
for  this  problem  (and  in  fact  a  more  general  class  of  problems) 
with  relatively  little  effort. 

Before  describing  the  general  assumptions  on  which 
these  results  are  based  it  will  be  useful  to  define  for  the 
seismic  case  described  above  a  few  parameters  of  interest. 

p[sgn(s+n^)  =  +1]  =  p[s+ni>0]  =  1-q  (4.1) 

Therefore , 

00  2  / 

q  =  /  -i-  e-X  /2  dx  =  Q(s/o)  (4.2) 

s/a  /27 

In  the  seismic  case  of  interest  we  expect  the  signals  to 
resemble  a  portion  of  a  sine  wave.  Furthermore,  for  detection, 
the  peak  signal  is  of  interest.  For  consistency  with  other 
discussions,  we  shall  therefore  define  SNR  as  if  s  were  the 
peak  value  of  a  sine  wave,  i.e., 

Is2 

SNR  =  (4.3) 

o 
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The  following  table  presents  some  pairs  of  values  of  SNR  and  q 
that  may  be  useful  in  subsequent  discussions. 

SNR  =  3.3  dB  -0.8  dB  -4.5  dB  -8.6  dB 

q  =  0.02  0.1  0.2  0.3 

4.1.1  Abstract  Mathematical  Model 

The  results  presented  below  apply  to  a  more  general 
problem  than  the  seismic  case  described  above.  To  illustrate 
this  point,  the  only  assumptions  that  are  crucial  to' the  develop¬ 
ment  will  now  be  stated. 

Based  on  the  observation  of  N  independent  random  varia¬ 
bles,  y^ ,  it  is  desired  to  determine  whether  hypothesis 
["signal  present"]  or  hypothesis  Hq  ["signal  not  present"] 
applies.  The  y^'s  are  described  by 

p[y,=+l|H  ]  =  l-p[y.=-l|H  ]  =  l-q>l/2 

1  1  11  (4.4) 

p[y1=+i|H()]  =  p[y1=-i  |hq]  =  1/2 

Based  on  these  assumptions  alone,  it  is  easy  to  show  that 

N 

v  =  l  y± 

i  1 

is  a  sufficient  statistic  and  that  the  optimum  decision  rule 
under  any  of  several  assumptions  is  equivalent  to  comparing 
the  sum  to  a  threshold,  K,  and  choosing  if  and  only  if  the 
threshold  is  exceeded  or  equalled. 

It  is  easily  seen  that  the  seismic  problem  described 
above  conforms  to  this  model,  and  it  could  be  modified  in  several 
ways  and  still  conform.  For  example,  provided  the  y^  are  inde¬ 
pendent  and  identically  distributed,  they  would  still  conform 
to  the  model  even  if  they  were  not  obtained  by  hardlimiting  the 
result  of  adding  a  Gaussian  random  variable  to  a  constant.  For 
example,  considering  an  experiment  with  N  independent  tosses  of 
a  coin,  y.  could  be  defined  as  +1  if  a  head  occurs  on  the  i-th 
toss  and  -1  if  a  tail  occurs,  and  the  model  would  still  apply. 

In  this  case,  the  hypothesis  Hq  would  correspond  to  the  coin 
being  "fair." 
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4.2  pf  ,  pfd  AND  THEIR  ASYMPTOTIC  FORMS 

4,2.1  False  Alarms 


The  probability  of  a  false  alarm  when  the  detection 
rule  suggested  above  is  employed  is  given  by 


Pfa  =  p[ u-K | HQ ]  (4.5) 

Denoting  the  number  of  negative  values  of  y^  by  X,  we  have 

U  =  N-2  X  (4.6) 


and  therefore 


fa 


r  -v  <  H — K  |  -f,  -1 

—  p  L  ^  ~  2  ^  q  J 

,  „  (N-KJ/2 

=  <j)H  l  (?) 

*  1=0  1 


(4.7) 


where  it  is  understood  that  N-K  is  an  even  number.  (This 

assumption  is  for  convenience  only  and  in  no  way  represents  a 

N — K 

restriction.  If  N-K  is  odd,  L— p— J ,  the  largest  integer  less 
I  —  K 

than  — replaces  the  latter.)  Since  Equation  (4.7)  can  be 
cumbersome  to  evaluate,  especially  when  N  is  large,  we  shall 
employ  upper  and  lower  bounds  that  were  originally  developed 
in  another  context.  From  Equation  (7.8)  of  Fano  [4]  and 
Equation  (4,7)  above,  we  have 


N  ,  M 

(  )(-)N 
''N-K''  '  2 1 


Pfa  < 


N+K 


N 


_ (  )(-)N 

2K^N-K' ' 2 ' 


(4.8) 


In  order  to  obtain  some  indication  of  the  tightness  of  these 
bounds,  it  should  be  noted  that  they  only  differ  by  a  factor 
of  N+K/2K.  If  for  example,  we  had  K  =  N/2,  the  factor  is  1.5. 

A  convenient  form  for  these  bounds  can  be  obtained 
by  considering  their  asymptotic  behavior  as  N  increases  and 
the  threshold  ratio,  r,  defined  by 


r  =  K/N 

remains  constant.  Defining 

-log2Pfa 


E  =  lim 

1  d  J\J  +» 


N 


(4.9) 


(4.10) 
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and  following  Equation  7.3  of  Fano  [4]  we  have 


„  ,  ,  1-r  ,  1-r  ,  1+r  ,  1+r 

Efa  =  1  +  log2-2—  +  ~2~  log- 
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(4.11) 


and 


-NE 


pfa  *  2 


fa 


(4.12) 


A  specific  example  may  help  to  illustrate  the  tight¬ 
ness  of  these  bounds.  Assume  N=21  and  K=17.  Then,  from 
Equation  (4.8) 

>-4 


1.1*10  4  <  p 


fa 


<  1.23*10 


Substituting  r  =  17/21  into  Equation  (4.11)  yields 
Efa  =  0.546 

and  therefore  from  Equation  (4.12) 

-4 


Pfa 


3*10 


(4.13) 


(4.14) 


(4.15) 


which  is  within  a  factor  of  3  of  the  correct  answer,  p~ 

_  ii  1  a 

1.11-10  . 


4.2.2 


False  Dismissals 


The  probability  of  a  false  dismissal  when  the  signal 
is  present  is  given  by 

Pfd  =  pCu<K  | 

r,  N-K  |  u 
=  pU>— —  |  H 1 J 


N 


i  (^)o-q 


N-i  i 


(4.16) 


i-S^+i 

where,  as  before,  A  is  the  number  of  negative  y^'s  and  N-K 
is  assumed  to  be  divisible  by  2.  Returning  again  to  Fano  [4], 
Equation  7.18,  we  have 


/  1+r. 

q(-o-) 


( l-q) (~2~  +  £) 


/  1+Ts 

q(-p-) 

B  <  Pfd  <  i_r" —  B  for  r<1_2q  (4.17) 


-q 
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where 


N+K  N-K 


N 


B  =  (N-K)(1-q)  2  q  2 


(4.18) 


As  In  the  case  of  false  alarms,  we  may  define 

log2pfd 


E„,  =  lim  - 
f  d  .. 

N->“ 


N 


(4.19) 


which,  following  Equations  7.44  and  7-58  of  Fano  [4],  may  be 
written  as 


Efd  2  l0S2  2q  +  2  l0g2  2( 1-q)  f0r  r<1  2q  ^*20) 

Finally,  we  have 

P_NEfd  (4.21) 

pfd 


4.2.3  Tradeoffs  Between  pfa  and  pfd 

It  Is  customary  In  detection  problems  to  Illustrate 
possible  tradoffs  between  p^a  and  p^d  by  plotting  receiver 
operating  characteristic  (ROC)  curves.  In  these  plots  the 
detection  probability  (1-p^)  Is  shown  as  a  function  of  Pfa> 
with  the  SNR  as  a  parameter.  For  large  N,  the  exact  expressions 
for  p^a  and  p...  (Equations  4.7  and  4.l6)  are  cumbersome  to 
calculate,  and,  furthermore,  all  of  the  calculations  must  be 
repeated  If  a  new  value  of  N  Is  considered.  By  employing 
the  exponents  E^a  and  E^,  Introduced  above,  most  of  this 
effort  can  be  avoided,  and  satisfactory  estimates  still  obtained. 

In  Figure  4.1  a  family  of  ROC-like  curves  for  E^d  and 
E^a  are  presented.  Each  curve  Is  obtained  by  varying  the 
threshold  parameter,  r,  while  holding  q  fixed  at  the  Indicated 
value.  (The  SNR's  that  correspond  to  these  q's  In  the  case 
of  signal  plus  Gaussian  noise  appear  In  the  table  Immediately 
preceding  Section  4.1.1.) 
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FIGURE  4.1 

DETECTION  CHARACTERISTICS 
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ASYMPTOTIC 


4.3 


COMPARISON  OF  DIMUS  AND  ANALOG  ARRAYS  IN  ADDITIVE 
GAUSSIAN  NOISE 


4.3.1  Properties  of  Analog  Array 


The  effects  of  hardlimiting  in  the  case  of  additive 
noise  have  been  considered  in  a  number  of  papers  [see,  for 
example,  references  5,6].  The  emphasis  of  these  papers  has 
usually  been  on  comparing  the  output  SNR's  of  the  DIMUS  and 
the  corresponding  analog  array.  Unfortunately,  the  results 
of  these  comparisons  are  not,  in  general,  the  answer  to  the 
following  question:  how  much  additional  SNR  at  the  input 

of  the  DIMUS  array  is  necessary  to  achieve  the  same  p^a  and 
p^  as  the  corresponding  analog  array?  It  is  this  question 
that  we  wish  to  answer. 

To  consider  this  question  we  first  need  expressions 
for  p^a  and  p^d  in  the  analog  case.  These  are,  in  fact,  easier 
to  obtain  than  in  the  DIMUS  case.  Making  only  a  slight  modi¬ 
fication  of  previous  definitions,  we  assume  that  the  array 
output,  u ,  is  compared  to  a  threshold  K 0.  Then 


fa 


N 

=  PC  I 

i=l 


n. 


L  Ko 


Hn]  =  Q(— ) 
/N 


(4.22) 


Upper  and  lower  bounds  on  p^a  can  be  obtained  by  using  the 
following  bounds  on  Q(*)  given  by  Equation  2.121  of  Reference  [3]. 


e"T  /2(1-  <  Q(y)  <  — —  e  y  72  y>0  (4.23) 

/2  TT  Y  Y  /2  TT  Y 

If,  as  before,  we  hold  r  =  K/N  fixed  and  let  N  increase,  then 
the  argument  of  Q(*)  in  Equations  ( 4 . 22 ) - ( 4 . 2 3)  will  be  equal 
to  r/N,  the  upper  and  lower  bounds  will  have  the  same  exponent, 

Ef a 5  given  by 

log?Q( r/N )  „ 

Efa  =  lim - ^ -  =  |r  log2e  (4.24) 


As  before,  then 

-NEfa  -N[|r2log2e] 
pfa  >  2  =  2 


(4.25) 
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The  probability  of  a  false  dismissal  is  given  by 


N  N-  -K 

Pfd  =  p[  l  (s+n1)<Ko]  =  — )  (4.26) 

i  =  l  /N 

In  the  same  manner  that  Equation  (4.24)  was  obtained,  we  may 
obtain 

Efd  =  ^(f  ~r)2  log2e,  s/a  >r  (4.27) 

Combining  Equations  (4.24)  and  (4.27)  yields 

Efd  =  (f/(log2e)/2  -  /E^)2  (4.28) 

The  asymptotic  form  of  p^  may  now  be  written  as 

-N  L  A  s?2oO  (log„e)-  /E“] 2 

Pfd  -  2  2  la  (4.29) 


4.3.2 


Comparison  of  DIMUS  and  Analog  Arrays 


Three  different  ways  of  viewing  the  difference  between 
the  DIMUS  and  analog  arrays  in  terms  of  the  exponents  are 
presented  in  Figure  4.2-4. 4.  Figure  4.2  presents  the  ROC-like 
curves  introduced  earlier  for  both  arrays  and  for  SNR's  of 
-0.8  dB  and  -8.6  dB.  Figure  4.3  considers  a  fixed  E^a  and  plots 
SNR  vs.  E^.  Finally,  Figure  4.4  plots  the  additional  SNR 
necessary  for  the  DIMUS  array  to  achieve  the  same  E^  as  the 
analog  array,  as  a  function  of  r,  or,  equivalently,  as  a 
function  of  E^a .  Strictly  speaking.  Figure  4.4  is  based  on  the 
limiting  case  of  E^d-»-0,  and  this  point  merits  a  few  additional 
comments.  Returning  to  Figure  4.3,  it  is  clear  that  the  SNR 
difference  between  the  DIMUS  and  analog  arrays  is  a  function 
of  E^^.  In  general  applications,  however,  we  are  principally 
interested  in  small  values  of  For  example,  with  N  =  100 

elements,  a  pfd  of  0.1  implies  an  E^d  of  0.03;  whereas  a  p„a 
of  10“6  would  imply  an  E^a  of  0.20.  As  suggested  by  this 
example  we  expect  to  be  concerned  with  values  of  E^  near  zero, 
and  values  of  E^a  that  are  significantly  larger.  With  this  in 
mind,  we  shall  take  as  a  standard  of  comparison  the  SNR  difference 
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FIGURE  4.2 

COMPARISON  OF  DETECTION  CHARACTERISTICS 
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FIGURE  4.3 

COMPARISON  OF  SNR'S  FOR 
E  f  a  AND  Efd 
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as  (Figure  4.3  suggests  that  this  will  give  a  smaller 

difference  than  applies  at  any  nonzero  value  of  E^.) 

The  Ef  =  0.278  point  of  Figure  4.4  is  obtained  by 
measuring  the  difference  (2.6  dB)  between  the  curves  of 
Figure  4.3  at  E^  =  0.  The  other  points  could,  in  principle, 
be  obtained  similarly,  although  in  fact  it  is  easier  to 
develop  equations  specifically  for  this  purpose.  From 
Figure  4.4  it  is  clear  that  the  SNR  loss  of  DIMUS  over  the 
analog  case  increases  with  K  from  a  minimum  value  of  just 
under  2  dB ,  which  is  exactly  the  tt/2  factor  in  power  that  is 
often  quoted  as  the  loss  due  to  hardlimiting . 

Figures 4 . 2-4 . 4  have  the  advantage  of  providing 
relationships  that  are  independent  of  N;  they  have  the 
disadvantage,  however,  that  some  additional  calculations 
are  necessary  to  go  from  the  exponents  (assuming  a  particular 
N)  to  the  more  familiar  quantities  pfa  and  pf^. 

Figures  4.5  and  4.6  present  pd  =  l-pfd  vs  •  SNR  curves 
for  several  assumed  values  of  N,  K  and  Pfa,  for  both  the 
analog  and  DIMUS  arrays.  These  curves  indicate  in  a  more 
concrete  fashion  some  of  the  conclusions  hinted  at  earlier 
and  illustrate  typical  values  that  result  when  parameters 
that  are  thought  to  be  relevant  to  seismic  applications  are 
employed.  It  should  be  noted  that  the  motives  for  the 
choices  of  N  =  21,  80  and  127  are  somewhat  mixed.  Some 
experiments  have  previously  been  performed  using  the  21 
center  elements  of  LASA-Montana.  Additional  experiments  are 
being  planned  that  may  use  approximately  80  elements  of 
LASA-Montana,  chosen  to  minimize  noise  correlation  between 
elements.  Finally,  we  understand  that  LASA-Norway  will 
eventually  include  in  excess  of  127  short  period  instruments. 
N=127  is  the  largest  number  that  could  be  quickly  and  easily 
handled  with  the  computer  and  subroutines  available  for  these 
illustrative  calculations. 

The  number" 2  dB"is  often  quoted  as  the  degradation 
resulting  from  bandlimiting .  Figure  4.4  suggests  that  the 
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PROBABILITY  OF  DETECTION 


FIGURE  4.5 

PROBABILITY  OF  DETECTION  CURVES 
(  P fa  ~  I0“3) 
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number  applies  to  a  limited  range  of  cases  if  probabilities 
are  considered.  Figure  4.6  makes  this  point  explicit.  The 
DIMUS  and  analog  curves  are  presented  for  the  N=21  case. 

For  DIMUS  the  false  alarm  probability  is  10  ,  for  the  analog 

array  it  is  somewhat  lower.  If  the  analog  threshold  were 
adjusted  to  yield  the  same  p^a,  the  curve  would  move  left  on 
the  figure.  Even  in  its  present  position,  however,  it  is 
clear  that  there  is  more  than  a  2  dB  discrepancy  between  the 
curves,  especially  at  the  higher  detection  probabilities. 

4.4  RELATIONSHIP  OF  DETECTION  RESULTS  TO  CODING  THEORY  RESULTS 

It  is  no  accident  that  the  bounds  provided  in  Refer¬ 
ence  [4]  are  useful  in  simplifying  the  results  of  Section  4.2. 

Those  familiar  with  bounds  on  the  probability  of  error  for  a 
binary  symmetric  channel  (BSC)  will,  after  some  interpretation, 
recognize  that  p^d  is  equal  to  the  probability  of  error  given 
by  the  "sphere  packed"  bound.  This  fundamental  bound  of  infor¬ 
mation  theory  provides  a  lower  bound  on  the  probability  of 
error  for  the  binary  symmetric  channel.  However,  for  the 
detection  problem,  the  sphere  packed  bound  is  not  a  bound  at 
all,  but  is  exactly  equal  to  the  probability  of  a  false  dismissal. 

A  list  of  the  relationships  between  the  detection  and 
coding  parameters  is  given  in  Table  V.  It  is  significant  to  note 
that  the  trade-off  in  coding  theory  between  the  value  of  the 
sphere  packed  exponent  and  the  rate  of  communication,  is  exactly 
analogous  to  the  trade-off  between  the  values  of  E^d  and  E^a 
given  by  Figure  4.1.  The  point  where  E^a  achieves  its  maximum 
value,  as  E^+0,  is  of  significance  in  coding  theory  since  it 
is  equal  to  the  channel  capacity  of  a  BSC.  It  is  also  of 
significance  in  detection  theory  since  it  corresponds  to  the 
trade-off  which  is  of  most  practical  interest. 

The  analog  between  detection  theory  and  coding  theory 
does  not  stop  with  the  relationship  between  an  array  of  hard- 
limited  data  and  data  passed  through  a  BSC.  Indeed,  a  similar 
relationship  exists  between  the  bounds  for  the  array  of  analog 
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TABLE  V 


RELATIONSHIP  BETWEEN  DETECTION  AND 


Parameter 


Detection  Theory 


N  The  number  of  elements  in 

the  array. 

q  The  probability  of  re¬ 

ceiving  a  -1. 

The  probability  of  a 
false  dismissal. 


The  false  dismissal 
exponent . 


max 


'fd 


E 


fa 


The  maximum  false  dis¬ 
missal  exponent  occurring 


when 


Efa>0 


The  false  alarm  exponent 


max  E  „ 
fa 


The  maximum  false  alarm 
exponent  occurring  when 


Efd"°- 


CODING  PARAMETERS 
Coding  Theory 


The  number  of  letters 
in  each  code  word 

The  crossover  proba¬ 
bility  for  the  BSC. 

The  probability  of  an 
error  for  the  sphere 
packed  bound 

The  sphere  packing 
error  exponent . 

The  zero  rate  sphere 
packing  experiment 


The  rate  of  communi¬ 
cation  . 

The  channel  capacity 
of  a  binary  symmetric 
channel . 
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data  and  the  results  known  about  the  Infinite  bandwidth  Gaussian 
channel.  It  seems  evident  that  coding  theory  techniques  developed 
for  more  complex  communications  channels,  those  involving  K>2 
symbols,  can  be  exploited  in  the  treatment  of  more  complex  detection 
problems,  and  perhaps,  the  reverse  is  also  true. 


-58- 


SECTION  V 


AUTOMATIC  pP  TEST 

An  automatic  pP  test  that  employs  seismograms  obtained 
from  a  continental-size  seismic  array  has  been  reported  in  pre¬ 
vious  reports  [1,2,7].  Recently,  the  properties  of  this  test 
have  been  re-examined  with  the  objective  of  answering  more 
clearly  two  questions:  under  what  circumstances  does  this 
test  perform  better  than  other  methods  and  what  components  of 
the  test  are  essential  to  its  success?  To  consider  these 
questions,  we  have  reconsidered  previously  processed  data  from 
a  variety  of  seismic  events:  earthquakes  with  depths  ranging 
from  5  km  to  150  km  and  underground  nuclear  explosions. 

The  automatic  test  is  composed  of  three  basic  parts: 

(1)  a  procedure  for  aligning  seismograms  for  coherent  addition 
of  the  pP  phase  arrivals;  (2)  an  energy  measurement  for  detec¬ 
tion  of  the  pP  phase  arrival  on  the  phased  sum  seismogram;  and 
(3)  a  measure  of  signal  correlation  among  the  seismograms  to 
distinguish  the  pP  phase  arrival,  which  would  be  fairly  consis¬ 
tent  across  the  array,  from  spurious  high  energy  peaks  on 
individual  records.  The  test  has  been  described  in  more  detail 
elsewhere  [2,7].  A  critical  examination  will  be  given  these 
components  of  the  automatic  test  in  an  attempt  to  justify  their 
necessity.  If  any  is  not  required  for  good  test  performance, 
then  it  may  be  possible  to  simplify  the  automatic  test  signi¬ 
ficantly.  This  is  the  case  for  the  shallow  pP  test  which 
searches  for  event  depths  in  the  10  to  40  km  region.  In  this 
region,  it  is  only  necessary  to  make  a  gross  correlation  measure 
as  defined  by  the  coda-correlation  discriminant  [Section  IX  of 
Reference  17],  and  then,  if  this  correlation  measure  indicates 
that  the  event  is  in  the  10  to  40  km  region,  to  calculate  simple 
energy  ratios  without  correlation  weighting.  Such  a  simplifica¬ 
tion  is  not  evident  for  the  extended  test  which  searches  for 
event  depths  down  to  200  km.  The  value  of  pP  phasing  will  be 
demonstrated  with  three  moderate  depth  earthquakes  (about  50  km). 
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The  correlation  measure  alone  has  proved  not  to  be  a  suitable 
statistic  for  detection  of  pP,  but  its  value  for  reducing 
"false"  energy  peaks  will  be  illustrated. 

Consideration  will  also  be  given  to  the  need  for 
widely  spaced  (continental-size)  seismometer  arrays,  rather 
than  moderate-size  arrays,  such  as  LASA-Montana.  The  value  of 
the  larger  arrays  for  reducing  local  noise,  aftershocks,  and 
local-to-the-seismometer  crustal  reverberations  in  the  pP 
phased  seismogram  and  in  calculations  of  the  n  test  statistic, 
will  be  noted. 

Finally,  in  previous  applications  of  the  test,  spurious 
peaks  which  tend  to  make  confident  depth  estimates  difficult, 
if  not  impossible,  occasionally  occur.  For  these  events 
erroneous  conclusions  may  result  as  is  the  case  of  a  nuclear 
explosion  in  the  Novaya  Zemlya  region.  These  extraneous  peaks 
can  be  explained  by  a  multilayered  model  of  the  earth's  crust. 
While  we  know  of  no  way  to  eliminate  this  problem,  it  does 
appear  that  a  knowledge  of  the  geology  of  the  epicenter  region 
may  be  very  useful  in  interpreting  the  test  results. 

5.1  ARRAY  PROCESSING  FOR  pP  ENHANCEMENT 

It  is  well-known  that  array  processing  can  be  used  to 
enhance  desired  signals  in  the  presence  of  undesirable  "noise" 
when  dealing  with  signals  more  strongly  correlated  than  the 
surrounding  noise.  This  is  generally  the  case  with  waveforms 
of  teleseismic  events.  The  pP  phase  arrival  is  the  signal  of 
interest,  and  the  surrounding  energy  is  considered  "noise." 
Portions  of  the  coda,  however,  can  be  as  strongly  correlated 
as  the  pP  signal  across  an  array. 

Earthquake  aftershocks,  crustal  reverberations  and 
local  noise  often  make  identification  of  the  pP  phase  arrival 
difficult  on  individual  seismograms.  An  array  such  as  LASA- 
Montana  can  alleviate  the  local  noise  problem,  but  local 
crustal  reverberations  in  the  coda  can  be  quite  correlated. 

This  is  evidenced  by  experimental  runs  of  the  automatic  pP  test 
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with  LASA  data.  When  the  individual  seismograms'  signal-to- 
noise  ratios  were  weak,  the  case  where  the  automatic  test  is  of 
most  interest,  multiple  peaks  occurred  in  the  test  statistic 
[Section  V  of  Reference  1],  probably  due  to  crustal  reverberations 
beneath  the  seismometers.  Larger,  continental-size  arrays  can 
eliminate  local  (to  the  seismometer)  crustal  reverberations 
since  the  earth's  crustal  structure  is  different  below  each 
receiver.  At  the  same  time,  these  larger  arrays  help  eliminate 
aftershocks.  With  large  arrays,  the  relative  P-pP  delay  times 
can  vary  by  as  much  as  6  seconds  for  a  200  km  depth  event,  whereas 
aftershocks  follow  main  P  by  a  constant  delay  time  at  all  seis¬ 
mometers,  provided  the  aftershock  hypocenter  is  the  same  as  the 
initial  shock  hypocenter.  Hence,  as  signals  are  phased  for  pP 
by  delaying  them  relative  to  P,  the  aftershock  signals  will  not 
add  coherently.  Crustal  reverberations  local  to  the  event  are 
another  problem,  and  will  be  discussed  in  a  later  section. 

To  demonstrate  the  power  of  array  processing  for  enhance¬ 
ment  of  the  pP  phase  arrival,  and  hence  to  show  the  advantage 
the  automatic  test  has  over  techniques  which  utilize  only  indi¬ 
vidual  seismograms,  it  was  felt  that  clear  visual  evidence  was 
required.  With  this  goal  in  mind,  a  technique  was  developed 
to  combine  a  set  of  seismograms  for  one  event  from  widely  spaced 
seismometers  into  one  phased-sum  seismogram  in  such  a  fashion 
that,  regardless  of  the  event's  depth  of  focus,  the  pP  phase 
arrivals  on  the  individual  records  would  add  coherently  in  the 
sum  seismogram.  For  a  given  depth  of  focus,  the  P-pP  delay 
time  increases  nonlinearly  with  range.  For  a  fixed  range  with 
a  constant  velocity  assumption  the  P-pP  delay  time  increases 
linearly  with  depth.  In  forming  the  sum  seismogram  the  seis¬ 
mometer  nearest  the  event  is  taken  as  a  standard.  Then  the 
seismograms  from  the  more  distant  stations  are  compressed  in 
time  by  deleting  points  so  that  the  P-pP  delay  times  are  com¬ 
parable  to  the  P-pP  delay  time  at  the  nearest  station.  Since 
the  delay  times  increase  linearly  with  depth,  points  are  removed 
from  seismograms  at  equally  spaced  intervals  and  the  remaining 
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points  shifted  together  as  if  the  points  removed  were  never 
there.  For  typical  arrays,  this  means  that  at  most  about  one 
point  in  ten  is  removed.  Hence,  the  compression  should  intro¬ 
duce  negligible  waveform  distortion  on  the  majority  of  the 
seismograms.  The  "compressed"  seismograms  are  then  simply 
summed  by  aligning  on  main  P,  and  the  pP  phases  should  add 
coherently,  regardless  of  the  true  depth  of  the  event. 

Three  events  are  presented  where  the  pP  (or  sP)  phase 
cannot  be  located  easily  on  individual  seismograms,  but  can  be 
located  on  the  phased  sum  seismogram.  The  first  event,  6  October 
1962,  was  an  Ecuador  earthquake  reported  earlier  [Section  VI 
of  Reference  2]  as  difficult  to  analyze;  our  tentative  results 
indicated  a  77  km  depth.  C&GS,  however,  had  very  good  evidence 
from  nearby  S-phase  arrivals  that  the  event  was  about  150  km 
deep.  Figure  5.1  shows  the  individual  seismograms  for  this 
event,  which  have  been  filtered  with  a  0.6-3.25  Hz  bandpass 
filter.  It  is  not  possible  to  identify  with  any  confidence 
the  pP  or  sP  phase  arrival  on  any  of  these  seismograms  even 
with  a  priori  information  that  this  is  a  150  km  event.  Figure 
5 . 3A  shows  the  phased  sum  of  the  "compressed"  versions  of  the 
seismograms  shown  in  Figure  5.1.  Here,  there  is  clearly  a  phase 
arrival  about  41.5  seconds  following  the  P-phase  arrival.  The 
nearest  station  is  4l.6°  from  the  event;  this  corresponds  well 
to  the  P-sP  delay  time  for  a  150  km  earthquake.  It  is  not 
unusual  for  the  sP  phase  to  be  much  stronger  than  the  pP  phase, 
which  seems  to  be  the  case  for  this  event. 

The  second  event,  2  January  1964,  was  a  magnitude-5.1 
Kamchatka  earthquake,  reported,  when  only  the  shallow  pP  test 
existed,  as  deeper  than  40  km  [Section  X  of  Reference  7].  The 
unfiltered  seismograms  for  this  event  are  shown  in  Figure  5.2. 

On  three  of  the  five  records,  the  pP  phase  arrival  is  difficult 
to  identify,  or  its  delay  time  relative  to  P  is  difficult  to 
measure.  In  addition  to  strong  noise,  there  are  strong  codas 
on  these  records.  Notice  what  happens,  however,  on  the  phased 
sum  of  the  "compressed"  seismograms  shown  in  Figure  5-3B.  The 
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FIGURE  5.1  -  SEISMOGRAMS  FOR  6  OCTOBER  1962  EVENT  (FILTERED) 
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FIGURE  5.2- SEISMOGRAMS  FOR  2  JANUARY  1964  EVENT 
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FIGURE  5.3  SUM  SEISMOGRAMS 


pP  phase  arrival  is  clear  and  the  general  coda  level  is  well 
down  relative  to  that  of  the  individual  seismograms.  It  is 
easy  to  make  a  confident  depth  estimate  using  this  phased  sum 
record.  The  measured  P-pP  delay  time  is  11.4  seconds;  the 
nearest  station  range  is  52.4°.  This  gives  a  depth  estimate 
of  42  km,  which  compares  well  with  the  40  km  estimate  of  the 
Vela  Seismological  Center.  Discussion  of  a  third  example, 
the  21  March  1962  Alaska  earthquake  is  deferred  until  the 
next  section,  since  it  is  used  to  illustrate  several  additional 
points,  and  fits  more  logically  into  that  discussion. 

The  two  events  just  discussed  illustrate  the  pP  or 
sP  enhancement  power  of  array  processing  designed  to  suppress 
interfering  coda  and  noise  which  might  be  strong  enough  to 
make  identification  on  individual  seismograms  difficult. 

A  reasonable  question  at  this  point  is  whether  or  not  these 
phase  arrivals  would  show  up  just  as  clearly  on  the  sum  of 
"uncompressed"  seismograms.  The  next  section  considers  this 
question. 

5.2  ALIGNMENT  OF  SEISMOGRAMS  FOR  THE  pP  PHASE  ARRIVAL 

In  order  to  utilize  the  coherent  property  of  the.  pP 
phase  arrivals  at  the  elements  of  a  seismometer  array,  the 
individual  signals  must  be  added  in  phase.  This  means  the 
relative  P-pP  delay  times  vs.  depth  must  be  known,  or  calcu¬ 
lated,  to  within  about  0.2  seconds  for  each  pair  of  array  ele¬ 
ments.  As  was  pointed  out  in  an  earlier  report  [Section  VI  of 
Reference  2],  this  level  of  precision  appears  feasible  with 
the  delay  time  calculations  used  by  the  automatic  test.  There 
are  indications,  however,  that  severe  anomalies  in  P-pP  delay 
times  are  not  unusual.  The  Seismic  Data  Laboratory  reports 
up  to  0.5  second  variation  between  the  P-pP  delay  time  indicated 
by  JB  tables  and  that  measured  at  the  TFO  Extended  Array  [8]. 
Since  the  TFO  Extended  Array  has  only  a  4°  aperture  and  the 
automatic  test  used  arrays  with  10°  to  30°  apertures,  delay 
time  discrepancies  could  be  even  more  severe  than  0.5  second. 
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This  seems  not  to  be  the  case  though,  at  least  on  a  majority  of 
the  seismograms  used  in  the  following  discussion. 

Experiments  indicate  that  P-pP  delay  times  can  be 
calculated  precisely  enough  to  enhance  significantly  the 
appearance  of  the  pP  phase  arrival  on  a  sum  seismogram  phased 
for  pP  as  compared  to  a  sum  seismogram  formed  simply  by  adding 
seismograms  phased  on  main  P.  Three  events  are  presented  to 
substantiate  this  result. 

The  sum  seismogram  without  pP  phasing  for  the  first 
event,  the  2  January  1964  event  of  the  last  section,  is  shown 
in  Figure  5 . 3C .  The  pP  phase  arrival  is  indicated  on  the  plot. 
This  should  be  compared  to  Figure  5.3B,  which  is  the  sum 
seismogram  phased  for  pP .  The  increased  clarity  of  Figure  5 • 3B 
is  evident. 

The  sum  seismogram  without  pP  phasing  for  the  second 
event,  the  15  July  1963  Kamchatka  earthquake,  is  shown  in 
Figure  5- 3D.  Below  this  is  the  sum  seismogram  phased  for  pP, 
Figure  5  •  3E .  Again,  the  pP  phase  arrival  is  signficantly  more 
well-defined  on  the  sum  seismogram  that  is  phased  for  pP. 

The  third  event  presented  illustrates  several  points. 
This  is  an  Alaskan  earthquake  of  21  March  1962.  Figure  5-4A 
shows  the  individual  siesmograms  for  the  event.  Figure  5-4B 
shows  the  sum  seismogram  with  phasing  for  pP.  Figure  5.4c 
shows  the  sum  seismogram  without  phasing  for  pP.  The  strong 
signal  following  main  P  by  55  seconds  appears  more  clearly 
in  the  sum  seismogram  without  pP  phasing  than  in  the  sum  seismo¬ 
gram  with  pP  phasing.  It  seems  this  is  an  aftershock.  Both 
we  and  C&GS  have  reported  this  "aftershock"  as  the  pP  phase 
arrival  [Section  VI  of  Reference  2].  The  relative  delays 
used  to  form  the  pP  phased  seismogram  have  decorrelated  the 
aftershock  waveforms  on  the  individual  records,  hence  the  mud¬ 
dled  appearance  of  this  phase  arrival  in  Figure  5.4B.  Simul¬ 
taneously,  the  pP  phasing  has  enhanced  the  true  pP  phase  arrival 
as  indicated  in  Figure  5.4B,  while  this  arrival  is  not  clear 
in  Figure  5.4c.  This  is  an  example  of  how  a  large  array  (1100  km 
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aperture)  can  alleviate  the  aftershock  problem.  The  n  statistic 
was  larger  for  the  pP  phase  arrival  than  for  the  aftershock 
arrival,  even  though  the  aftershock  clearly  shows  more  energy 
than  the  pP  phase.  This  event  also  illustrates  the  point  of  the 
first  section:  when  the  pP  phase  cannot  be  identified  clearly 
on  individual  seismograms,  it  may  be  distinguishable  on  the  sum 
seismogram.  It  should  be  noted  that  this  event  uses  only  four 
stations,  a  marginal  number  for  good  automatic  test  performance 
with  the  strong  coda  of  the  records  used.  These  four  stations, 
however,  have  a  good  "average"  spread  in  range,  that  is,  all 
stations  are  separated  by  at  least  300  km.  Also,  the  array 
is  reasonably  close  to  the  event,  the  station  ranges  varying 
from  about  25°  to  36°.  In  this  region,  the  travel-time  tables 
have  their  greatest  slope,  giving  maximum  relative  P-pP  delay 
times  for  an  array  with  a  given  aperture. 

5-3  ENERGY  AND  CORRELATION  DETECTION  OP  pP 

The  development  of  the  automatic  pP  test  was  an  ad  hoc 
process.  The  test  began  by  measuring  the  energy  in  a  one- 
second  window  at  the  assumed  location  of  the  pP  phase  arrival 
on  the  sum  seismogram  and  comparing  this  to  the  surrounding 
average  energy  level.  The  test  was  applied  only  to  the  10  to 
^0  km  depth  region,  and  peak  values  of  the  calculated  energy 
ratio  were  selected  as  the  pP  phase  arrivals.  It  was  reasoned 
at  that  time  that  spurious  high  energy  peaks  on  one  seismogram 
could  possibly  change  the  location  of  the  energy  ratio  peak 
from  the  true  pP  location  since  the  summing  technique  tended 
to  emphasize  seismograms  with  high  signal-to-noise  ratios.  Thus 
the  average  paired  correlation  of  the  signals  in  the  trial  loca¬ 
tions  of  the  pP  arrivals  on  the  individual  seismograms  was  intro 
duced  as  a  weighting  factor,  with  the  hope  of  eliminating  these 
spurious  energy  peaks.  In  retrospect,  it  appears  this  corre¬ 
lation  measure  was  not  required  for  this  purpose.  In  2k  earth¬ 
quake  trial  runs,  the  correlation  weighting  factor  did  not  once 
change  the  location  of  the  test  statistic  peak,  that  is,  the 
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energy  ratio  and  the  energy  ratio  weighted  by  the  average 
paired  correlation  peaked  at  the  same  place  for  each  trial  run. 
However,  in  comparing  energy  ratio  peaks  of  earthquakes  with 
those  of  nuclear  explosions,  one  finds  roughly  equivalent  results. 
If  a  reasonable  energy  ratio  threshold  is  determined  for  the  24 
earthquake  runs  in  such  a  fashion  that  peaks  of  the  energy  ratio 
above  this  threshold  are  accepted  as  valid  pP  phase  arrivals, 
and  peaks  below  the  threshold  value  are  rejected  as  spurious 
energy  peaks,  then  for  100%  detection  of  valid  pP  arrivals,  14 
out  of  22  runs  which  should  not  indicate  a  pP  phase  arrival 
in  the  10  to  40  km  region  would  have  energy  ratios  large  enough 
to  exceed  this  threshold;  11  of  these  14  "false"  energy  peaks 
occur  with  shots.  If  now  the  average  paired  correlation  is 
introduced  to  weight  the  energy  ratio  and  a  new  threshold  similar 
to  the  above  is  determined  for  the  weighted  energy  statistic, 
there  are  only  three  "false"  peaks,  two  of  them  from  shots, 
which  would  be  accepted  as  valid  pP  phase  arrivals.  The  one 
false  earthquake  peak  was  for  the  15  July  1963  Kamchatka  event, 
which  will  be  considered  along  with  the  two  shots,  in  the  dis¬ 
cussion  of  crustal  reverberations. 

A  coda  correlation  discriminant  was  developed  simulta¬ 
neously  with  the  shallow  version  (10  to  40  km)  of  the  automatic 
pP  test.  This  discriminant  is  quite  effective  at  separating 
10  to  40  km  earthquakes  from  other  events.  One  major  simpli¬ 
fication  of  the  automatic  pP  test  in  this  shallow  region  might 
be  to  run  first  the  coda-correlation  test  to  determine  if  an 
event  is  to  be  accepted  as  "in  the  10  to  40  km  region."  If  so, 
then  the  energy  ratio  would  be  calculated  without  correlation 
weighting  and  the  peak  would  be  selected  as  the  pP  location. 

If  the  event  is  not  "in  the  10  to  40  km  region,"  according  to 
the  coda  correlation  test,  then  no  pP  test  need  be  run.  With 
this  scheme,  only  one  shot  and  one  earthquake  produced  "false" 
results . 
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Since  these  results  were  obtained,  the  pP  test  has  been 
extended  to  search  for  possible  event  depths  down  to  200  km. 

A  measure  similar  to  the  original  coda  correlation  discriminant 
has  not  yet  been  developed  for  deeper  events,  thus  the  pP  test 
still  has  some  significant  value  for  discrimination  purposes. 

For  deeper  test  runs,  it  is  possible  for  aftershocks  to  yield 
large  values  for  the  energy  ratio  statistic,  such  as  in  the 
Alaskan  earthquake  of  21  March  1962.  Also,  shots  may  yield 
spurious  energy  peaks  which  probably  could  be  eliminated  by 
correlation  weighting.  At  the  present  time,  not  enough  shots 
have  been  run  with  the  extended  test  to  determine  whether  or  not 
a  problem  of  spurious  peaks  in  the  energy  ratio  exists. 

There  has  been  no  previous  discussion  of  pP  detection 
using  only  a  correlation  measure  among  the  individual  seismograms. 
Based  on  previous  test  calculations,  it  is  possible  to  state  that 
such  a  test  statistic  would  not  give  satisfactory  results.  In 
37  experimental  runs,  using  correlation  alone  to  detect  pP, 
misleading  results  occur  in  13  cases.  "Misleading  results"  mean 
deviations  of  the  pP  locations  indicated  by  correlation  from 
the  presently  accepted  pP  locations  based  on  the  best  available 
information.  Of  15  shot  runs,  five  runs  gave  high  enough  corre¬ 
lation  values  to  be  mistaken  for  the  pP  phase  arrival.  Unex¬ 
pectedly  high  correlation  values  may  well  be  due  to  signal  re¬ 
flections  from  velocity  discontinuities  in  the  earth's  crust. 

Since  the  velocity  discontinuities,  for  example  8.0  km/sec  to 
6.3^  km/sec,  are  not  as  severe  as  the  earth-air  or  earth-water 
interface,  their  reflected  signals  will  not  contain  as  much 
energy,  in  general,  as  the  pP  or  sP  phases.  Hence  the  need  for 
an  energy  detection  scheme. 
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5.4 


THE  EFFECT  OF  VELOCITY  DISCONTINUITIES  IN  THE  EARTH'S  CRUST 


As  discussed  earlier,  several  interfering  types  of  signals 
can  be  adequately  handled  by  the  automatic  pP  test  by  using 
continental-size  arrays  with  seismometer  separations  of  at  least 
200  km.  Crustal  reverberations  local  to  the  event,  however,  are 
more  difficult  to  handle.  Earthquake  signal  reflections  will 
occur  at  velocity  discontinuities  within  the  earth's  crust,  such 
as  the  Mohorovicic  discontinuity.  Evidence  exists  to  show  the 
presence  of  other  discontinuities,  such  as  the  Conrad  discon¬ 
tinuity  occurring  at  a  depth  of  about  10  km  [9-13].  In  some  regions 
of  the  earth,  this  discontinuity  may  be  more  severe  in  terms  of 
the  ratio  of  velocities  on  either  side  than  the  Mohorovicic 
discontinuity.  Reflections  from  these  discontinuities  can 
produce  signals  which  are  coherent,  even  across  a  continental- 
size  array,  and  which,  furthermore,  have  delay  times  relative 
to  main  P  that  are  approximately  the  same  as  the  P-pP  delay 
times.  Four  experimental  runs  indicate  these  reflections  can 
significantly  affect  the  results  of  the  pP  test. 

The  first  event  affected  by  a  reflection  at  a  crustal- 
velocity  discontinuity  is  the  Kamchatka  earthquake  of  15  July 
1963.  This  event  was  reported  by  the  Vela  Seismological  Center 
(VSC)  as  a  60  km  earthquake.  When  the  shallow  pP  test  was  run, 
results  indicated  the  event  was  a  19  km  earthquake.  The  coda 
correlation  discriminant  showed  that  this  event  is  "in  the  10 
to  40  km  region."  [Appendix  A  of  Reference  7-]  On  applying 
the  extended  automatic  pP  test,  results  then  confirmed  the  VSC 
depth  estimate  of  60  km,  which  was  further  verified  by  inspec¬ 
tion  of  the  original  seismograms.  This  "false"  peak  at  19  km 
is  probably  due  to  the  earthquake  signal  reflection  from  the 
Moho,  which  seems  to  be  about  35  km  deep  in  this  region.  The 
"true"  peak  at  60  km  exhibits  an  energy  ratio  of  12  compared  to 
an  energy  ratio  of  2.4  at  19  km.  The  phases  reflected  from 
the  Moho  are  a  bit  more  correlated  though,  enough  so  as  to 
fool  the  coda-correlation  test. 
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The  second  event  affected  by  such  discontinuities  is 
the  Sea  of  Okhotsk  earthquake  of  25  August  1963.  In  the  shallow 
test  results,  three  strong  peaks  occurred  in  the  n  statistic, 
one  at  12  km,  a  second  at  22  km,  and  a  third  at  30  km.  This 
could  be  accounted  for  by  a  30  km  event,  with  reflections  from 
the  Moho  and  Conrad  discontinuities.  A  detailed  picture  of  the 
crustal  structure  in  the  Okhotsk  region  is  given  by  Healy  [12], 
and  the  calculated  delay  times  for  such  reflections  are  comparable 
with  results  of  the  automatic  pP  test. 

The  last  two  events  affected  are  shots,  one  in  the 
Novaya  Zemlya  region  and  one  in  the  Kazakh  region.  Both  events 
gave  n  peaks  of  greater  than  1.0  when  the  shallow  pP  test  was 
applied.  The  Novaya  Zemlya  peak  occurred  at  a  test  depth  of 
15  km,  which  could  have  been  the  result  of  a  signal  path  from 
the  event  down  to  the  Conrad  discontinuity  (at  an  assumed  depth 
of  10  km)  back  to  the  surface  and  on  to  the  seismometer.  The 
Kazakh  peak  occurred  slightly  deeper,  about  18  km,  which  could 
have  been  caused  by  similar  reflections  with  a  somewhat  deeper 
Conrad  discontinuity  or  a  somewhat  slower  propagation  velocity. 

The  Novaya  Zemlya  reflections  were  so  correlated  as  to  fool  the 
coda  correlation  discriminant  test.  The  Kazakh  event  did  not 
do  so . 

Problems  with  these  interfering  reflection  arrivals 
are  inherent  to  the  test  as  it  is  presently  designed.  Crustal 
reflection  delay  times  should  be  at  least  roughly  comparable 
to  P-pP  delay  times,  and  since  these  signals  come  so  soon  after 
the  P  arrival,  the  relative  P-pP  delay  times  will  not  have  a 
great  variation  across  even  a  continental-size  array.  If  these 
reflections  show  up  consistently  and  strongly  with  shots,  such 
as  they  may  do  in  the  Novaya  Zemlya  region,  it  may  be  necessary 
to  study  crustal  structure  in  problem  areas  in  order  to  learn 
what  to  expect  if  a  nuclear  explosion  is  detonated. 
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5-5 


SUMMARY  AND  CONCLUSIONS 


A  retrospective  examination  of  the  automatic  pP  test 
has  been  made  in  order  to  ascertain  the  potential  utility  of  the 
test  to  make  reasonable  depth  estimates  when  more  standard 
techniques  based  on  inspection  of  individual  seismograms  may 
fail,  to  understand  possible  limitations  on  the  test's  perfor¬ 
mance,  and  to  consider  possible  simplifications  in  the  test's 
implementation.  A  few  of  the  more  difficult  events  which  have  been 
run  with  the  automatic  test  have  been  examined  in  some  detail. 

The  power  of  array  processing  designed  for  pP  enhance¬ 
ment  has  been  verified  to  some  extent  with  three  examples,  the 
6  October  1962  Ecuador  earthquake,  the  2  January  1964  Kamchatka 
earthquake,  and  the  21  March  1962  Alaska  earthquake.  For  all 
three  events,  either  identification  of  the  pP  (or  sP)  phase 
arrivals,  or  measurement  of  their  delay  times  relative  to  P 
was  difficult  on  individual  seismograms,  yet  the  pP  (or  sP)  phase 
arrival  was  fairly  easily  distinguished  on  a  sum  seismogram 
formed  to  emphasize  the  pP  (or  sP)  phase  arrival.  Use  of 
continental-size  arrays  for  pP  processing  eliminates  the  problem 
of  local-to-the-seismometers  crustal  reverberations. 

To  apply  array  processing  for  pP  enhancement  it  must  be 
possible  to  add  the  individual  signals  in  phase,  which  means 
relative  P-pP  delay  times  must  be  predictable  for  a  given  range 
and  depth  to  within  about  0.2  second.  Although  another  report 
indicates  that  discrepancies  of  measured  P-pP  times  from  J-B 
predicted  values  by  as  much  as  0.5  second  may  not  be  unusual, 
three  examples  have  been  presented  which  demonstrate  that  phasing 
of  seismometer  signals  for  pP  enhancement  offers  a  significant 
improvement  in  detection  capability  over  seismometer  signals 
summed  by  simply  phasing  on  main  P.  Using  continental-size 
arrays,  these  phasing  times  may  be  as  large  as  6  seconds  for  a 
200  km  event.  Such  large  delay  times  can  eliminate  the  problem 
of  mistaking  aftershocks  for  pP. 

By  reviewing  previous  test  results,  it  has  been  found 
that  the  shallow  pP  test  (10  to  40  km  depths)  performs  as  well 
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with  just  an  energy  ratio  as  it  does  with  an  energy  ratio  weighted 
by  a  correlation  measure  for  detection  of  pP.  However,  energy 
ratios  are  comparable  for  10-40  km  earthquakes  and  shots.  Corre¬ 
lation  weighting  allows  the  dismissal  of  "false"  energy  peaks 
for  all  but  three  of  the  twenty-two  runs  on  which  no  peak  should 
appear  in  the  10-40  km  region.  If  instead  of  correlation  weighting 
of  the  energy  ratio  statistic,  a  gross  correlation  measure  (based 
on  a  10  sec  time  window)  is  calculated,  a  simplified  pP  test  could 
be  used  and  only  two  classification  errors  in  twenty-two  runs  would 
result.  A  similar  simplification  for  the  extended  automatic  pP 
test  is  not  evident  at  this  time.  Thus  the  present  design  of  the 
extended  test,  that  is  aligning  for  pP,  and  then  calculating  both 
energy  ratios  and  correlation  measures  will  be  retained  in  the 
near  future,  at  least  until  more  experimental  runs  are  performed 
with  shots. 

Velocity  dis continuities  in  the  earth's  crust  were  shown 
to  have  affected  the  automatic  test  results  on  four  runs,  three 
of  which  were  misclassified  by  the  shallow  pP  test  and  two  of 
which  were  misclassified  by  the  coda  correlation  discriminant . 

The  problems  are  believed  to  have  resulted  from  signal  reflections 
from  either  the  Mohorovicic  discontinuity,  generally  about  30  km 
deep,  or  the  Conrad  dis continuity ,  10  km  deep.  Signals  reflected 
from  these  dis continuities  can  contain  significant  energy  and 
will  be  correlated  even  across  continental  arrays.  They  can  also 
have  delay  times  relative  to  main  P  comparable  with  P-pP  delay 
times.  So  far  this  has  presented  an  uncorrectable  problem  for 
the  automatic  test.  One  approach  to  resolving  matters  may  be  to 
study  crustal  structure  in  problem  areas  to  learn  what  to  expect 
when  an  event  occurs.  This  may  permit  proper  interpretation  of 
test  results. 
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APPENDIX  A 


WEAK  NOISE  ERROR  IN  CORRELATION  TIME  PICK 


We  assume  a  signal  at  time  t'  of  the  form 

C  (  i.  4-  I  ^  —  ,/p  OU  ®  2  Tt  W  (  t-t  '  ) 

S  1 5 1  )  -  /ES2W  2^  (t-t') 


( A-l ) 


where  W  is  the  bandwidth  of  the  signal. 

If  the  value  of  t'  is  to  be  estimated,  the  situation 
resembles  the  communication  problem  of  pulse  position  modula¬ 
tion  (PPM) . 

Assuming  the  noise  to  be  Gaussian  with  flat  spectrum 
Nq/2  in  the  same  band  as  the  signal,  we  have  from  Reference  [3], 


2  V2 
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where 

S2  =  /  Ijfr  S( t ,t ’ ) | 2  dt 
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Thus 
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This  integral  is  better  evaluated  in  the  frequency 
domain,  e.g.,  by  Parseval's  theorem. 

The  spectrum  of  i|>(t)  is 


♦  (f) 
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elsewhere 
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Differentiation  in  time  corresponds  to  multiplication  by  j2irf 
in  frequency,  thus 
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and  therefore, 
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